基于并行EBE-CG方法的有限元求解程序（从我的毕业论文中节选出来的）

附录A  EBE-CG法主程序
/* FemPar.c
  Written by XiaoFei,college of mathematics science and computing technology,CSU,2009-5  */
#include<stdio.h>
#include<malloc.h>
#include"mpi.h"
#define PROCNO 4
typedef struct
{
    float Pos[2];
    int Local;
    int Type;
}NodeType;
typedef struct
{
    int Vertex[3];
}ElementType;
NodeType *Node;
ElementType *Element;
int *Shared,MaxCommon,Common[PROCNO],*Neighbours[PROCNO];
int ProcNo,ProcID,Nodes,Elements,IntNodes,IBNodes;
float *Buf;
 
float InnerProduct(float *A1,float *A2,float *B1,float *B2)
{
    int I;
    float IP=0.0,IPTocal;
    for(I=0;I<IntNodes;I++)
        IP+=A1[I]*B1[I];
    for(I=0;I<IBNodes;I++)
        IP+=(A2[I]*B2[I])/Shared[I];
    MPI_Allreduce(&IP,&IPTocal,1,MPI_FLOAT,MPI_SUM,MPI_COMM_WORLD);
    return(IPTocal);
}
void Update(float *A)
{
    int I,J;
    MPI_Status Status;
    for(I=0;I<PROCNO;I++)
    {
        for(J=0;J<Common[I];J++)
            Buf[J]=A[Node[Neighbours[I][J]].Local];
        if(Common[I]>0)
            MPI_Send(&Buf[0],Common[I],MPI_FLOAT,I,10,MPI_COMM_WORLD);
    }
    for(I=0;I<PROCNO;I++)
    {
        if(Common[I]>0)
            MPI_Recv(&Buf[0],Common[I],MPI_FLOAT,I,10,MPI_COMM_WORLD,&Status);
        for(J=0;J<Common[I];J++)
            A[Node[Neighbours[I][J]].Local]+=Buf[J];
    }
}
int CG(float **Ap,float **As,float **Bp,float *Fp,float *Fs,float *Vp,float *Vs)
{
    float *Rp,*Rs,*Pp,*Ps,*Qp,*Qs,GammaOld,GammaNew,Tau,Alpha,Beta,Tol=1.0e-4;
    int I,J,K,KMax=250;
    Pp=(float*)malloc(IntNodes*sizeof(float));
    Ps=(float*)malloc(IBNodes*sizeof(float));
    Qp=(float*)malloc(IntNodes*sizeof(float));
    Qs=(float*)malloc(IBNodes*sizeof(float));
    Rp=(float*)malloc(IntNodes*sizeof(float));
    Rs=(float*)malloc(IBNodes*sizeof(float));
    Update(Fs);
    for(I=0;I<IntNodes;I++)
        Pp[I]=Rp[I]=Fp[I];
    for(I=0;I<IBNodes;I++)
        Ps[I]=Rs[I]=Fs[I];
    GammaNew=InnerProduct(Rp,Rs,Rp,Rs);
    if(GammaNew<Tol*Tol)
        return(1);
    if(ProcID==0)
        printf("Gamma = %f/n",GammaNew);
    for(K=0;K<KMax;K++)
    {
        for(I=0;I<IntNodes;I++)
        {
            for(J=0,Qp[I]=0.0;J<IntNodes;J++)
                Qp[I]+=Ap[I][J]*Pp[J];
            for(J=0;J<IBNodes;J++)
                Qp[I]+=Bp[I][J]*Ps[J];
        }
        for(I=0;I<IBNodes;I++)
        {
            for(J=0,Qs[I]=0.0;J<IntNodes;J++)
                Qs[I]+=Bp[I][J]*Pp[J];
            for(J=0;J<IBNodes;J++)
                Qs[I]+=As[I][J]*Ps[J];
        }
        Update(Qs);
        Tau=InnerProduct(Pp,Ps,Qp,Qs);
        Alpha=GammaNew/Tau;
        for(I=0;I<IBNodes;I++)
        {
            Vp[I]+=Alpha*Pp[I];
            Rp[I]-=Alpha*Qp[I];
        }
        for(I=0;I<IBNodes;I++)
        {
            Vs[I]+=Alpha*Ps[I];
            Rs[I]-=Alpha*Qs[I];
        }
        GammaOld=GammaNew;
        GammaNew=InnerProduct(Rp,Rs,Rp,Rs);
        if(ProcID==0)
            printf("Gamma = %f (K = %d)/n",GammaNew,K);
        if(GammaNew<Tol*Tol)
            return(1);
        Beta=GammaNew/GammaOld;
        for(I=0;I<IntNodes;I++)
            Pp[I]=Rp[I]+Beta*Pp[I];
        for(I=0;I<IBNodes;I++)
            Ps[I]=Rs[I]+Beta*Ps[I];
    }
    return(0);
}
float BC(float X,float Y)
{
    return(X+Y);
}
 
main(int argc,char** argv)
{
    FILE *InFile,*OutFile;
    char FileName[40];
    int I,J,K,II,JJ;
    float **Ap,**Bp,**As,*Fp,*Fs,*Vp,*Vs;
    float SLoc[3][2],TwoA,Area,Gr0[3],Gr1[3],StiffJI;
 
    MPI_Init(&argc,&argv);
    MPI_Comm_rank(MPI_COMM_WORLD,&ProcID);
    MPI_Comm_size(MPI_COMM_WORLD,&ProcNo);
    if(ProcNo!=PROCNO)
    {
        printf("Unexpected number of processors./n");
        MPI_Abort(MPI_COMM_WORLD,-1);
    }
    if(ProcID<10)
        sprintf(FileName,"Data0%d.In",ProcID);
    else
        sprintf(FileName,"Data%d.In",ProcID);
    InFile=fopen(FileName,"r");
    fscanf(InFile,"%d%d",&Nodes,&Elements);
   
    Node=(NodeType*)malloc(Nodes*sizeof(NodeType));
    for(I=0,IntNodes=0,IBNodes=0;I<Nodes;I++)
    {
        fscanf(InFile,"%f%f%d",&Node[I].Pos[0],&Node[I].Pos[1],&Node[I].Type);
        if(Node[I].Type==1)
            Node[I].Local=IntNodes++;
        else if(Node[I].Type>1)
           Node[I].Local=IBNodes++;
        else
           Node[I].Local=-1;
    }
    Shared=(int*)malloc(IBNodes*sizeof(int));
    for(I=0,IBNodes=0;I<Nodes;I++)
        if(Node[I].Type>1)
            Shared[IBNodes++]=Node[I].Type;
    Element=(ElementType*)malloc(Elements*sizeof(ElementType));
    for(I=0;I<Elements;I++)
        fscanf(InFile,"%d%d%d",&Element[I].Vertex[0],&Element[I].Vertex[1],&Element[I].Vertex[2]);
    MaxCommon=0;
    for(I=0;I<PROCNO;I++)
    {
        fscanf(InFile,"%d",&Common[I]);
        if(Common[I]>MaxCommon)
            MaxCommon=Common[I];
        if(Common[I]>0)
            Neighbours[I]=(int*)malloc(Common[I]*sizeof(int));
        for(J=0;J<Common[I];J++)
            fscanf(InFile,"%d",&Neighbours[I][J]);
    }
    Buf=(float*)malloc(MaxCommon*sizeof(float));
    fclose(InFile);
    printf("ProcID = %d,IntNodes = %d,IBNodes = %d/n",ProcID,IntNodes,IBNodes);
    Ap=(float**)calloc(IntNodes,sizeof(float*));
    for(I=0;I<IntNodes;I++)
        Ap[I]=(float*)calloc(IntNodes,sizeof(float));
    Bp=(float**)calloc(IBNodes,sizeof(float*));
    for(I=0;I<IntNodes;I++)
        Bp[I]=(float*)calloc(IntNodes,sizeof(float));
    As=(float**)calloc(IBNodes,sizeof(float*));
    for(I=0;I<IBNodes;I++)
        As[I]=(float*)calloc(IBNodes,sizeof(float));
    Fp=(float*)calloc(IntNodes,sizeof(float));
    Fs=(float*)calloc(IBNodes,sizeof(float));
    Vp=(float*)calloc(IntNodes,sizeof(float));
    Vs=(float*)calloc(IBNodes,sizeof(float));
    for(K=0;K<Elements;K++)
    {
        SLoc[0][0]=Node[Element[K].Vertex[0]].Pos[0];
        SLoc[1][0]=Node[Element[K].Vertex[1]].Pos[0];
        SLoc[2][0]=Node[Element[K].Vertex[2]].Pos[0];
        SLoc[0][1]=Node[Element[K].Vertex[0]].Pos[1];
        SLoc[1][1]=Node[Element[K].Vertex[1]].Pos[1];
        SLoc[2][1]=Node[Element[K].Vertex[2]].Pos[1];
        TwoA=SLoc[0][0]*(SLoc[1][1]-SLoc[2][1])+
             SLoc[1][0]*(SLoc[2][1]-SLoc[0][1])+
             SLoc[2][0]*(SLoc[0][1]-SLoc[1][1]);
        Area=0.5*TwoA;
        Gr0[0]=(SLoc[1][1]-SLoc[2][1])/TwoA;
        Gr0[1]=(SLoc[2][1]-SLoc[0][1])/TwoA;
        Gr0[2]=(SLoc[0][1]-SLoc[1][1])/TwoA;
        Gr0[0]=-(SLoc[1][0]-SLoc[2][0])/TwoA;
        Gr0[1]=-(SLoc[2][0]-SLoc[0][0])/TwoA;
        Gr0[2]=-(SLoc[0][0]-SLoc[1][0])/TwoA;
        for(J=0;J<3;J++)
        {
            JJ=Element[K].Vertex[J];
            if(Node[JJ].Type==1)
                for(I=0;I<3;I++)
                {
                    II=Element[K].Vertex[I];
                    StiffJI=Area*(Gr0[I]*Gr0[J]+Gr1[I]*Gr1[J]);
                    if(Node[II].Type==1)
                        Ap[Node[JJ].Local][Node[II].Local]+=StiffJI;
                    else
                        if(Node[II].Type>1)
                            Bp[Node[JJ].Local][Node[II].Local]+=StiffJI;
                        else
                            if(Node[II].Type==0)
                            Fp[Node[JJ].Local]-=BC(Node[II].Pos[0],Node[II].Pos[1])*StiffJI;
                }
            else
                if(Node[JJ].Type>1)
                    for(I=0;I<3;I++)
                    {
                        II=Element[K].Vertex[I];
                        StiffJI=Area*(Gr0[I]*Gr0[J]+Gr1[I]*Gr1[J]);
                        if(Node[II].Type>1)
                            As[Node[JJ].Local][Node[II].Local]+=StiffJI;
                        else
                            if(Node[JJ].Type==0)
                            Fs[Node[JJ].Local]-=BC(Node[II].Pos[0],Node[II].Pos[1])*StiffJI;
                    }
        }
    }
    CG(Ap,As,Bp,Fp,Fs,Vp,Vs);
    for(I=0;I<Nodes;I++)
        if(Node[I].Type==1)
            printf(" %f %f %f/n",Node[I].Pos[0],Node[I].Pos[1],Vp[Node[I].Local]);
        else
            if(Node[I].Type>1)
                printf(" %f %f %f/n",Node[I].Pos[0],Node[I].Pos[1],Vs[Node[I].Local]);
    MPI_Finalize();
}
 
 
  
 
 
 
 
附录B  任务分配预处理程序
#include<stdio.h>
#define FILENAME "mesh.txt"
#define PROCNO 4
typedef struct{
float Pos[2]; //位置坐标
int GlobalDOF;
int LocalNode;
int Type;
}NodeType;
typedef struct{
    int Vertex[3];
}ElementType;
NodeType* Node;
ElementType* Element;
int PartSize[PROCNO],*PartE[PROCNO],*PartN[PROCNO];
 
main(int argc,char** argv)
{
    FILE *InFile,*OutFile;
    int Nodes,Elements,I,J,K,V,NodeCounter[PROCNO],Common[PROCNO];
    char FileName[40];
   
    InFile=fopen(FILENAME,"r");
//读入节点信息
    fscanf(InFile,"%d",&Nodes);//节点数目
    Node=(NodeType*)malloc(Nodes*sizeof(NodeType));
    for(I=0;I<PROCNO;I++)
        PartN[I]=(int*)malloc(Nodes*sizeof(int));
    for(I=0;I<Nodes;I++)
    {
      fscanf(InFile,"%d%f%f",&Node[I].GlobalDOF,&Node[I].Pos[0],&Node[I].Pos[1]);
        //读入节点GlobalDOF，坐标，
        Node[I].Type=0;
        //初始化节点符号，<0表示边界
    }
//读入单元信息
    fscanf(InFile,"%d",&Elements); //单元数目
    Element=(ElementType*)malloc(Elements*sizeof(ElementType));
    for(I=0;I<PROCNO;I++)
    {
        PartE[I]=(int*)malloc(Elements*sizeof(int));
        PartSize[I]=0;
    }
    for(I=0;I<Elements;I++)
        fscanf(InFile,"%d%d%d",&Element[I].Vertex[0],&Element[I].Vertex[1],&Element[I].Vertex[2]);
        //读入单元的三个vertex
//构造PartE
    for(J=0;J<Elements;J++)
    {
        fscanf(InFile,"%d",&I);
        PartE[I][PartSize[I]++]=J;
    }
    fclose(InFile);
//关闭文件
    for(I=0;I<PROCNO;I++)
    {
        NodeCounter[I]=0;
        for(J=0;J<Nodes;J++)
        {
            Node[J].LocalNode=-1;
            PartN[I][J]=0;
        }//初始化
//构造PartN
        for(J=0;J<PartSize[I];J++)
        {
            for(K=0;K<3;K++)
            {
                V=Element[PartE[I][J]].Vertex[K];
                if(Node[V].LocalNode==-1)
                    Node[V].LocalNode=NodeCounter[I]++;
                if(Node[V].GlobalDOF>-1&&PartN[I][V]==0)
                {
                    Node[V].Type++;
                    PartN[I][V]=1;
                }
            }
        }
    }
    for(I=0;I<PROCNO;I++)
    {
        if(I<10)
            sprintf(FileName,"Data0%d.In",I);
        else
            sprintf(FileName,"Data%d.In",I);
        OutFile=fopen(FileName,"w");
        fprintf(OutFile," %d %d/n",NodeCounter[I],PartSize[I]);
        NodeCounter[I]=0;
        for(J=0;J<Nodes;J++)
            Node[J].LocalNode=-1;
        for(J=0;J<PartSize[I];J++)
        {
            for(K=0;K<3;K++)
            {
                V=Element[PartE[I][J]].Vertex[K];
                if(Node[V].LocalNode==-1)
                {
                    Node[V].LocalNode=NodeCounter[I]++;
                    fprintf(OutFile," %f %f %d/n",Node[V].Pos[0],Node[V].Pos[1],Node[V].Type);
                }
            }
        }
        for(J=0;J<PartSize[I];J++)
        {
            for(K=0;K<3;K++)
            {
                V=Element[PartE[I][J]].Vertex[K];
                fprintf(OutFile," %d",Node[V].LocalNode);
            }
            fprintf(OutFile,"/n");
        }
        for(K=0;K<PROCNO;K++)
        {
            Common[K]=0;
            if(I!=K)
                for(J=0;J<Nodes;J++)
                    if(PartN[I][J]>0&&PartN[K][J]>0)
                        Common[K]++;
        }
        for(K=0;K<PROCNO;K++)
        {
            fprintf(OutFile," %d/n",Common[K]);
            if(Common[K]>0)
            {
                for(J=0;J<Nodes;J++)
                    if(PartN[I][J]>0&&PartN[K][J]>0)
                        fprintf(OutFile," %d",Node[J].LocalNode);
                    fprintf(OutFile,"/n");
            }
        }
        fclose(OutFile);
    }
}