求解压力备份（）

2013-04-11之前的版本

求解压力poisson方程

#include <iomanip>
#include <fstream>
#include "matvec.h"
#include"GS.h"
#include"GMRES.h"
#include"ggje.h"
#include "callumfpack.h"
//2013-2-22 核查了A和B是否于Poisson离散出来形式一致
//加入了左边界条件为强制边界(Dirichlet条件)
//#define Temp(i,j) Temp[(i)*(N+1)*(M+1)+j] //使用宏定义将二维转换成一维
//2013-4-2 添加了umfpack计算压力
void Sol_P(int N,int M,double delta_x,double delta_y,double delta_t ,double **ustar ,double **vstar ,double **upie ,double **vpie,double **Pres)
{

    //void (*p)(double * ,int , double [],double []); //定义一个函数指针变量 p
    //cout<<"N="<<N<<" "<<"M="<<M<<endl;
    double **A=new double *[(N+1)*(M+1)];
    for(int k=0;k<(N+1)*(M+1);k++)
        A[k]= new double [(N+1)*(M+1)];
    for(int i=0;i<(N+1)*(M+1);i++)
        for(int j=0;j<(N+1)*(M+1);j++)
            A[i][j]=0;
    //cout<<endl;
    // 输出初始化的 A矩阵
    // cout<<" 初始化A矩阵 "<<endl;
    /*     for(int i=0;i<(N+1)*(M+1);i++)
    {
    for(int j=0;j<(N+1)*(M+1);j++)
    {
    cout<<A[i][j]<<" ";
    }
    cout<<endl;
    }*/
    double *B =new double [(N+1)*(M+1)];
    for(int bk=0;bk<(N+1)*(M+1);bk++)
    {
        B[bk]=0;
        //             cout<<B[bk]<<" ";
    }
    //cout<<endl;
    double *sol_P =new double [(N+1)*(M+1)];
    for(int bk=0;bk<(N+1)*(M+1);bk++)
    {
        sol_P[bk]=0;
        //     cout<<sol_P[bk]<<" ";
    }
    //cout<<endl;

    /*---------------------------------------------------------------------------------------------*/
    
    /*点(0,0)处的离散方程*/
    A[0*(M+1)+0][0*(M+1)+0]=-2*(1/pow(delta_x,2)+1/pow(delta_y,2));
    A[0*(M+1)+0][0*(M+1)+1]=2*(1/pow(delta_y,2));
    A[0*(M+1)+0][1*(M+1)+0]=2*(1/pow(delta_x,2));
    B[0*(M+1)+0]=1/delta_t*(ustar[0+1][0]-ustar[0][0])/delta_x+1/delta_t*(vstar[0][0+1]-vstar[0][0])/delta_y+2/delta_x*(-(upie[0][0]-ustar[0][0])/delta_t)+2/delta_y*(-(vpie[0][0]-vstar[0][0])/delta_t);
    //                                    向后差分
    /*i=0上的点（除去(0,0),(0,M)）*/
    for(int j=1;j<=M-1;j++)
    {
        int i=0;//第一行
        A[i*(M+1)+j][i*(M+1)+j-1]=1/pow(delta_y,2);
        A[i*(M+1)+j][i*(M+1)+j]=-2*(1/pow(delta_x,2)+1/pow(delta_y,2));
        A[i*(M+1)+j][i*(M+1)+j+1]=1/pow(delta_y,2);
        A[i*(M+1)+j][(i+1)*(M+1)+j]=2.0/pow(delta_x,2);
        B[i*(M+1)+j]=1/delta_t*(ustar[i+1][j]-ustar[i][j])/delta_x+1/delta_t*(vstar[i][j+1]-vstar[i][j-1])/(2*delta_y)+2/delta_x*(-(upie[i][j]-ustar[i][j])/delta_t);
        //                                   向后差分和中心差分（能用中心差分就用中心差分）
    }
    /*点(0,M)处的离散方程*/
    A[0*(M+1)+M][0*(M+1)+M-1]=2*(1/pow(delta_y,2));
    A[0*(M+1)+M][0*(M+1)+M]=-2*(1/pow(delta_x,2)+1/pow(delta_y,2));
    A[0*(M+1)+M][1*(M+1)+M]=2*(1/pow(delta_x,2));
    B[0*(M+1)+M]=1/delta_t*(ustar[0+1][M]-ustar[0][M])/delta_x+1/delta_t*(vstar[0][M]-vstar[0][M-1])/delta_y+2/delta_x*(-(upie[0][M]-ustar[0][M])/delta_t)-2/delta_y*(-(vpie[0][M]-vstar[0][M])/delta_t);
    /*i=1~N-1,j=1~M-1 之间的离散方程*/
    for(int i=1;i<N;i++)
    {
        for(int j=1;j<M;j++)
        {
            A[i*(M+1)+j][i*(M+1)+j-1]=1*(1/pow(delta_y,2));
            A[i*(M+1)+j][i*(M+1)+j]=-2*(1/pow(delta_x,2)+1/pow(delta_y,2));
            A[i*(M+1)+j][i*(M+1)+j+1]=1*(1/pow(delta_y,2));
            A[i*(M+1)+j][(i-1)*(M+1)+j]=1*(1/pow(delta_x,2));
            A[i*(M+1)+j][(i+1)*(M+1)+j]=1*(1/pow(delta_x,2));
            B[i*(M+1)+j]=1/delta_t*((ustar[i+1][j]-ustar[i-1][j])/(2*delta_x)+(vstar[i][j+1]-vstar[i][j-1])/(2*delta_y));//对一阶导数采用中心差分
        }
    }
    /*------------------------------------------------------------------------------------------------------*/
    /*点(N,0)处的离散方程*/
    A[N*(M+1)+0][N*(M+1)+0]=-2*(1/pow(delta_x,2)+1/pow(delta_y,2));
    A[N*(M+1)+0][N*(M+1)+1]=2*(1/pow(delta_y,2));
    A[N*(M+1)+0][(N-1)*(M+1)+0]=2*(1/pow(delta_x,2));
    B[N*(M+1)+0]=1/delta_t*(ustar[N][0]-ustar[N-1][0])/delta_x+1/delta_t*(vstar[N][0+1]-vstar[N][0])/delta_y-2/delta_x*(-(upie[N][0]-ustar[N][0])/delta_t)+2/delta_y*(-(vpie[N][0]-vstar[N][0])/delta_t);
    /*i=N上的点（除去(0,0),(0,M)）*/
    for(int j=1;j<=M-1;j++)
    {
        int i=N;//第N行
        A[i*(M+1)+j][i*(M+1)+j-1]=1/pow(delta_y,2);
        A[i*(M+1)+j][i*(M+1)+j]=-2*(1/pow(delta_x,2)+1/pow(delta_y,2));
        A[i*(M+1)+j][i*(M+1)+j+1]=1/pow(delta_y,2);
        A[i*(M+1)+j][(i-1)*(M+1)+j]=2.0/pow(delta_x,2);
        B[i*(M+1)+j]=1/delta_t*(ustar[i][j]-ustar[i-1][j])/delta_x+1/delta_t*(vstar[i][j+1]-vstar[i][j-1])/(2*delta_y)-2/delta_x*(-(upie[i][j]-ustar[i][j])/delta_t);
    }
    /*点(N,M)处的离散方程*/
    A[N*(M+1)+M][N*(M+1)+M-1]=2*(1/pow(delta_y,2));
    A[N*(M+1)+M][N*(M+1)+M]=-2*(1/pow(delta_x,2)+1/pow(delta_y,2));
    A[N*(M+1)+M][(N-1)*(M+1)+M]=2*(1/pow(delta_x,2));
    B[N*(M+1)+M]=1/delta_t*(ustar[N][M]-ustar[N-1][M])/delta_x+1/delta_t*(vstar[N][M]-vstar[N][M-1])/delta_y-2/delta_x*(-(upie[N][M]-ustar[N][M])/delta_t)-2/delta_y*(-(vpie[N][M]-vstar[N][M])/delta_t);
    //cout<<"B[N*(M+1)+M]="<<B[N*(M+1)+M]<<endl;
    /*---------------------------------------------------------------------------------------------------------------------*/
    /*j=M上的点（除去(0,M),(N,M)）*/
    for(int i=1;i<=N-1;i++)
    {
        int j=M;//第N行
        A[i*(M+1)+j][(i-1)*(M+1)+j]=1/pow(delta_y,2);
        A[i*(M+1)+j][i*(M+1)+j]=-2*(1/pow(delta_x,2)+1/pow(delta_y,2));
        A[i*(M+1)+j][i*(M+1)+j-1]=2*(1/pow(delta_y,2));
        A[i*(M+1)+j][(i+1)*(M+1)+j]=1.0/pow(delta_x,2);
        B[i*(M+1)+j]=1/delta_t*(ustar[i+1][j]-ustar[i-1][j])/(2*delta_x)+1/delta_t*(vstar[i][j]-vstar[i][j-1])/delta_y-2/delta_y*(-(vpie[i][j]-vstar[i][j])/delta_t);
    }
    /*---------------------------------------------------------------------------------------------------------------------*/
    /*j=0上的点（除去(0,0),(N,0)）*/
    for(int i=1;i<=N-1;i++)
    {
        int j=0;//第0行
        A[i*(M+1)+j][(i-1)*(M+1)+j]=1/pow(delta_x,2);
        A[i*(M+1)+j][i*(M+1)+j]=-2*(1/pow(delta_x,2)+1/pow(delta_y,2));
        A[i*(M+1)+j][i*(M+1)+j+1]=2*(1/pow(delta_y,2));
        A[i*(M+1)+j][(i+1)*(M+1)+j]=1.0/pow(delta_x,2);
        B[i*(M+1)+j]=1/delta_t*(ustar[i+1][j]-ustar[i-1][j])/(2*delta_x)+1/delta_t*(vstar[i][j+1]-vstar[i][j])/delta_y+2/delta_y*(-(vpie[i][j]-vstar[i][j])/delta_t);
    }
    /*---------------------------------------------------------------------------------------------------------------------*/
    //double *Temp=new double[((N+1)*(M+1))*((N+1)*(M+1))]; //动态分配内存空间用于存储临时矩阵 ;

    //将矩阵Temp全部元素初始化为
    //注：row为行号， col为列号
    
    
    /*-----------------------处理Ax=b中A的形式让其利用Gauss或者稀疏可解-----------------*/
    //????????????????????????????????
    for(int i=0;i<M+1;i++)
    {
        for(int j=0;j<(N+1)*(M+1);j++)
        {
            A[i][j]=0;
        }
        A[i][i]=1;
        B[i]=10;//左边界赋值为常数10
    }

    //ofstream outAMatrix;
    //outAMatrix.open("F:\\fluid\\C-V-IBM\\AMatrix.txt");
    ////速度太慢
    //for(int row=0;row<(N+1)*(M+1);row++)
    //{
    //    for(int col=0;col<(N+1)*(M+1);col++)
    //    {
    //        //Temp(row,col)=A[row][col];
    //        outAMatrix<<setprecision(5)<<A[row][col]<<" ";
    //    }
    //    outAMatrix<<endl;
    //}

    /*ofstream outBVector;
    outBVector.open("F:\\fluid\\C-V-IBM\\BVector.txt");
    for(int i=0;i<(N+1)*(M+1);i++)
    {
            outBVector<<setprecision(5)<<B[i]<<" ";
            outBVector<<endl;
    }*/
    double *testx =new double [(N+1)*(M+1)];
    for(int bk=0;bk<(N+1)*(M+1);bk++)
    {
        testx[bk]=0;
        //             cout<<testx[bk]<<" ";
    }
    //cout<<endl;
    /*---------------------调用UMFPACK求解 A*testx=B ---------------------------*/
    umf(A,B,testx,(N+1)*(M+1));
    for(int row=0;row<=N;row++)
    {
        for(int col=0;col<=M;col++)
        {
            //Pres[row][col]=sol_P[(row-1)*M+col-1];  //Pres represents Pressure
            Pres[row][col]=testx[row*M+col];// 调用umf  之前写成row*(M)+col 出现问题
        }
    }
    std::ofstream outPres;
    outPres.open("F:\\fluid\\C-V-IBM\\Pres.txt");    
    //cout<< "输出求解的矩阵Pres"<<endl;
    for(int i=0;i<N+1;i++)
    {
        for(int j=0;j<M+1;j++)
        {
            //cout<<Pres[i][j]<< " ";
            outPres<<setprecision(5)<<Pres[i][j]<<" ";
        }
        //cout<<endl;
        outPres<<endl;
    }
    
    for(int t=0;t<(N+1)*(M+1);t++)
        delete []A[t];
    delete []A;

    /*撤销u*空间*/

    delete [] B;
    delete [] sol_P;
    outPres.close();
    //outAMatrix.close();
}

这个有问题 主要有

1 压力不对称

2 压力输出有问题

Pres[row][col]=testx[row*M+col];的错误 困扰了很久

修改版本

#include <iomanip>
#include <fstream>
#include "matvec.h"
#include"GS.h"
#include"GMRES.h"
#include"ggje.h"
#include "callumfpack.h"
//2013-2-22 核查了A和B是否于Poisson离散出来形式一致
//加入了左边界条件为强制边界(Dirichlet条件)
//#define Temp(i,j) Temp[(i)*(N+1)*(M+1)+j] //使用宏定义将二维转换成一维
//2013-4-2 添加了umfpack计算压力
//2013-4-8 修改了在边界处离散的错误，以前为在每个边界按照一个方向差分 pian_u/pian_x=(u_n+1-u_n-1)/delta_x
//在边界处的差分都是由边界内部节点到外部节点
//2013-4-9 修改了pian_ustar/pian_x的迎风格式 按照具体物理模型进行离散
//2013-4-11 修改了压力输出的错误 以前导致有两条斜着的带子出现，错位导致
void Sol_P(int N,int M,double delta_x,double delta_y,double delta_t ,double **ustar ,double **vstar ,double **upie ,double **vpie,double **Pres)
{

    //void (*p)(double * ,int , double [],double []); //定义一个函数指针变量 p
    //cout<<"N="<<N<<" "<<"M="<<M<<endl;
    double **A=new double *[(N+1)*(M+1)];
    for(int k=0;k<(N+1)*(M+1);k++)
        A[k]= new double [(N+1)*(M+1)];
    for(int i=0;i<(N+1)*(M+1);i++)
        for(int j=0;j<(N+1)*(M+1);j++)
            A[i][j]=0;
    //cout<<endl;
    // 输出初始化的 A矩阵
    // cout<<" 初始化A矩阵 "<<endl;
    /*     for(int i=0;i<(N+1)*(M+1);i++)
    {
    for(int j=0;j<(N+1)*(M+1);j++)
    {
    cout<<A[i][j]<<" ";
    }
    cout<<endl;
    }*/
    double *B =new double [(N+1)*(M+1)];
    for(int bk=0;bk<(N+1)*(M+1);bk++)
    {
        B[bk]=0;
        //             cout<<B[bk]<<" ";
    }
    //cout<<endl;
    double *sol_P =new double [(N+1)*(M+1)];
    for(int bk=0;bk<(N+1)*(M+1);bk++)
    {
        sol_P[bk]=0;
        //     cout<<sol_P[bk]<<" ";
    }
    //cout<<endl;

    /*---------------------------------------------------------------------------------------------*/
    
    /*点(0,0)处的离散方程*/
    A[0*(M+1)+0][0*(M+1)+0]=-2*(1/pow(delta_x,2)+1/pow(delta_y,2));
    A[0*(M+1)+0][0*(M+1)+1]=2*(1/pow(delta_y,2));
    A[0*(M+1)+0][1*(M+1)+0]=2*(1/pow(delta_x,2));
    B[0*(M+1)+0]=1/delta_t*(ustar[0+1][0]-ustar[0][0])/delta_x+1/delta_t*(vstar[0][0+1]-vstar[0][0])/delta_y+2/delta_x*(-(upie[0][0]-ustar[0][0])/delta_t)+2/delta_y*(-(vpie[0][0]-vstar[0][0])/delta_t);
    //                                    向后差分
    /*i=0上的点（除去(0,0),(0,M)）*/
    for(int j=1;j<=M-1;j++)
    {
        int i=0;//第一行
        A[i*(M+1)+j][i*(M+1)+j-1]=1/pow(delta_y,2);
        A[i*(M+1)+j][i*(M+1)+j]=-2*(1/pow(delta_x,2)+1/pow(delta_y,2));
        A[i*(M+1)+j][i*(M+1)+j+1]=1/pow(delta_y,2);
        A[i*(M+1)+j][(i+1)*(M+1)+j]=2.0/pow(delta_x,2);
        B[i*(M+1)+j]=1/delta_t*(ustar[i+1][j]-ustar[i][j])/delta_x+1/delta_t*(vstar[i][j+1]-vstar[i][j-1])/(2*delta_y)+2/delta_x*(-(upie[i][j]-ustar[i][j])/delta_t);
        //                                   向后差分和中心差分（能用中心差分就用中心差分）
    }
    /*点(0,M)处的离散方程*/
    A[0*(M+1)+M][0*(M+1)+M-1]=2*(1/pow(delta_y,2));
    A[0*(M+1)+M][0*(M+1)+M]=-2*(1/pow(delta_x,2)+1/pow(delta_y,2));
    A[0*(M+1)+M][1*(M+1)+M]=2*(1/pow(delta_x,2));
    B[0*(M+1)+M]=1/delta_t*(ustar[0+1][M]-ustar[0][M])/delta_x+1/delta_t*(vstar[0][M]-vstar[0][M-1])/delta_y+2/delta_x*(-(upie[0][M]-ustar[0][M])/delta_t)-2/delta_y*(-(vpie[0][M]-vstar[0][M])/delta_t);
    /*i=1~N-1,j=1~M-1 之间的离散方程*/
    for(int i=1;i<N;i++)
    {
        for(int j=1;j<M;j++)
        {
            A[i*(M+1)+j][i*(M+1)+j-1]=1*(1/pow(delta_y,2));
            A[i*(M+1)+j][i*(M+1)+j]=-2*(1/pow(delta_x,2)+1/pow(delta_y,2));
            A[i*(M+1)+j][i*(M+1)+j+1]=1*(1/pow(delta_y,2));
            A[i*(M+1)+j][(i-1)*(M+1)+j]=1*(1/pow(delta_x,2));
            A[i*(M+1)+j][(i+1)*(M+1)+j]=1*(1/pow(delta_x,2));
            B[i*(M+1)+j]=1/delta_t*((ustar[i+1][j]-ustar[i-1][j])/(2*delta_x)+(vstar[i][j+1]-vstar[i][j-1])/(2*delta_y));//对一阶导数采用中心差分
        }
    }
    /*------------------------------------------------------------------------------------------------------*/
    /*点(N,0)处的离散方程*/
    A[N*(M+1)+0][N*(M+1)+0]=-2*(1/pow(delta_x,2)+1/pow(delta_y,2));
    A[N*(M+1)+0][N*(M+1)+1]=2*(1/pow(delta_y,2));
    A[N*(M+1)+0][(N-1)*(M+1)+0]=2*(1/pow(delta_x,2));
    B[N*(M+1)+0]=1/delta_t*(ustar[N][0]-ustar[N-1][0])/delta_x+1/delta_t*(vstar[N][0+1]-vstar[N][0])/delta_y-2/delta_x*(-(upie[N][0]-ustar[N][0])/delta_t)+2/delta_y*(-(vpie[N][0]-vstar[N][0])/delta_t);//改变下边界符号
    /*i=N上的点（除去(0,0),(0,M)）*/
    for(int j=1;j<=M-1;j++)
    {
        int i=N;//第N行
        A[i*(M+1)+j][i*(M+1)+j-1]=1/pow(delta_y,2);
        A[i*(M+1)+j][i*(M+1)+j]=-2*(1/pow(delta_x,2)+1/pow(delta_y,2));
        A[i*(M+1)+j][i*(M+1)+j+1]=1/pow(delta_y,2);
        A[i*(M+1)+j][(i-1)*(M+1)+j]=2.0/pow(delta_x,2);
        B[i*(M+1)+j]=1/delta_t*(ustar[i][j]-ustar[i-1][j])/delta_x+1/delta_t*(vstar[i][j+1]-vstar[i][j-1])/(2*delta_y)-2/delta_x*(-(upie[i][j]-ustar[i][j])/delta_t);
    }
    /*点(N,M)处的离散方程*/
    A[N*(M+1)+M][N*(M+1)+M-1]=2*(1/pow(delta_y,2));
    A[N*(M+1)+M][N*(M+1)+M]=-2*(1/pow(delta_x,2)+1/pow(delta_y,2));
    A[N*(M+1)+M][(N-1)*(M+1)+M]=2*(1/pow(delta_x,2));
    B[N*(M+1)+M]=1/delta_t*(ustar[N][M]-ustar[N-1][M])/delta_x+1/delta_t*(vstar[N][M]-vstar[N][M-1])/delta_y-2/delta_x*(-(upie[N][M]-ustar[N][M])/delta_t)-2/delta_y*(-(vpie[N][M]-vstar[N][M])/delta_t);
    //cout<<"B[N*(M+1)+M]="<<B[N*(M+1)+M]<<endl;
    /*---------------------------------------------------------------------------------------------------------------------*/
    /*j=M上的点（除去(0,M),(N,M)）*/
    for(int i=1;i<=N-1;i++)
    {
        int j=M;//第N行
        A[i*(M+1)+j][(i-1)*(M+1)+j]=1/pow(delta_y,2);
        A[i*(M+1)+j][i*(M+1)+j]=-2*(1/pow(delta_x,2)+1/pow(delta_y,2));
        A[i*(M+1)+j][i*(M+1)+j-1]=2*(1/pow(delta_y,2));
        A[i*(M+1)+j][(i+1)*(M+1)+j]=1.0/pow(delta_x,2);
        B[i*(M+1)+j]=1/delta_t*(ustar[i+1][j]-ustar[i-1][j])/(2*delta_x)+1/delta_t*(vstar[i][j]-vstar[i][j-1])/delta_y-2/delta_y*(-(vpie[i][j]-vstar[i][j])/delta_t);
    }
    /*---------------------------------------------------------------------------------------------------------------------*/
    /*j=0上的点（除去(0,0),(N,0)）*/
    for(int i=1;i<=N-1;i++)
    {
        int j=0;//第0行
        A[i*(M+1)+j][(i-1)*(M+1)+j]=1/pow(delta_x,2);
        A[i*(M+1)+j][i*(M+1)+j]=-2*(1/pow(delta_x,2)+1/pow(delta_y,2));
        A[i*(M+1)+j][i*(M+1)+j+1]=2*(1/pow(delta_y,2));
        A[i*(M+1)+j][(i+1)*(M+1)+j]=1.0/pow(delta_x,2);
        B[i*(M+1)+j]=1/delta_t*(ustar[i+1][j]-ustar[i-1][j])/(2*delta_x)+1/delta_t*(vstar[i][j+1]-vstar[i][j])/delta_y+2/delta_y*(-(vpie[i][j]-vstar[i][j])/delta_t);//改变下边界
    }
    /*---------------------------------------------------------------------------------------------------------------------*/
    //double *Temp=new double[((N+1)*(M+1))*((N+1)*(M+1))]; //动态分配内存空间用于存储临时矩阵 ;

    //将矩阵Temp全部元素初始化为
    //注：row为行号， col为列号
    //下面八行需要注意---需满足delta_x=delta_y
    for(int i=0;i<=N;i++)
    {
        for(int j=0;j<=M;j++)
        {
            A[i][j]=pow(delta_x,2.0)*A[i][j];
            B[i*(M+1)+j]=pow(delta_x,2.0)*B[i*(M+1)+j];
        }
    }
    
    /*-----------------------处理Ax=b中A的形式让其利用Gauss或者稀疏可解-----------------*/
    //????????????????????????????????
    for(int i=0;i<M+1;i++)
    {
        for(int j=0;j<(N+1)*(M+1);j++)
        {
            A[i][j]=0;
        }
        A[i][i]=1;
        B[i]=10;//左边界赋值为常数10
    }

    //ofstream outAMatrix;
    //outAMatrix.open("F:\\fluid\\cC-V-IBM\\AMatrix.txt");
    ////速度太慢
    //for(int row=0;row<(N+1)*(M+1);row++)
    //{
    //    for(int col=0;col<(N+1)*(M+1);col++)
    //    {
    //        //Temp(row,col)=A[row][col];
    //        outAMatrix<<setprecision(5)<<A[row][col]<<" ";
    //    }
    //    outAMatrix<<endl;
    //}

    ofstream outBVector;
    outBVector.open("F:\\fluid\\cC-V-IBM\\BVector.txt");
    /*for(int i=0;i<(N+1)*(M+1);i++)
    {
    outBVector<<setprecision(5)<<B[i]<<" ";
    outBVector<<endl;
    }*/
    for(int j=M;j>=0;j--)
    {
        for(int  i=0;i<N+1;i++)
        {
            //cout<<" "<<u[i][j]; 
            outBVector<<" "<<B[i*(M+1)+j];
            //outBVector<<endl;
        }
        //cout<<endl;
        outBVector<<endl;
    }
    outBVector.close();
    double *testx =new double [(N+1)*(M+1)];
    for(int bk=0;bk<(N+1)*(M+1);bk++)
    {
        testx[bk]=0;
        //             cout<<testx[bk]<<" ";
    }
    //cout<<endl;
    /*---------------------调用UMFPACK求解 A*testx=B ---------------------------*/
    umf(A,B,testx,(N+1)*(M+1));
    for(int row=0;row<=N;row++)
    {
        for(int col=0;col<=M;col++)
        {
            //Pres[row][col]=sol_P[(row-1)*M+col-1];  //Pres represents Pressure
            Pres[row][col]=testx[row*(M+1)+col];// 调用umf
        }
    }
    std::ofstream outPres;
    outPres.open("F:\\fluid\\cC-V-IBM\\Pres.txt");    
    //cout<< "输出求解的矩阵Pres"<<endl;
    for(int i=0;i<N+1;i++)
    {
        for(int j=0;j<M+1;j++)
        {
            //cout<<Pres[i][j]<< " ";
            outPres<<setprecision(5)<<Pres[i][j]<<" ";
        }
        //cout<<endl;
        outPres<<endl;
    }
    
    for(int t=0;t<(N+1)*(M+1);t++)
        delete []A[t];
    delete []A;

    /*撤销u*空间*/

    delete [] B;
    delete [] sol_P;
    outPres.close();
    //outAMatrix.close();
}