基于并行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);

    }

}

  • 0
    点赞
  • 2
    收藏
    觉得还不错? 一键收藏
  • 2
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论 2
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值