m.txt

c语言写得一阶有限元的通用程序

# include<stdio.h>
# include<stdlib.h>
# include<math.h>


# define LND 700
# define LEL 1200
#define PI 3.1415927
# define MU0 (4*PI*(1.0e-10))



struct Nodetype {float x;float y;};
struct elementtype {int i;int j;int m;};
struct boundaryNodetype {int bNode;float bptal;};


void KMatricStorageAddress();	/*Form the addressofK matric*/
void KPMatricSetup();			/* Form K & P matric*/
void BoundaryConditionDeal();	/* Dealwith boundary condition * /
void solveq();					/* Resolve thematric equation*/
void mbbub();					/* Sorting the coordinates*/


void solveq(a,m,b,id,n)
int m,n,*id;
float *a,*b;

{	int i,j,k,l,l1,l2,l3,k1,k2;

	for(i=0;i<n-1;i++)

	{l=id[i];

		for(j=i+1;j<n;j++)

		{	l1=j+1-(id[j]-id[j-1]);
			if(l1<=i)
			{	l2=id[j]-j+i;

				for(k=j;k<n;k++)
				{	l3=k+1-(id[k]-id[k-1]);

					if(l3<=i)
					{ 	k1=id[k]-k+i;
						k2=id[k]-k+j;
						a[k2]=a[k2]-a[k1]*a[l2]/a[l];  /* ???k1*/
					}
				}

				b[j]=b[j]-b[i]*a[l2]/a[l];

			}
		}
	}

	b[n-1]=b[n-1]/a[m -1];

	for(i=0;i<n-1;i++)

	{	k=n-2-i;

		for(j=k+1;j<n;j++)

		{	l1=j+1-(id[j]-id[j-1]);

			if(l1<=k)
			{l2=id[j]-j+k;b[k]=b[k]-b[j]*a[l2];
			}
		}

		l=id[k];
		b[k]=b[k]/a[l];
	}
}


main ()
{
	int k,nk,BeginNumber,SegmentTotal,LASTND,LASTEL,*Address;
	int BlockTotal,*BeNumber;
	float Rlarge,*Beta,*KMatric,*PMatric;
	struct Nodetype *Node;
	struct elementtype *Element;
	FILE *FilePoint,*fp;

	if((FilePoint=fopen("AGEN.dat","rb"))==NULL)

	{printf("cannotopen AGEN.datfile!\n");exit(0);}



	fread(&BlockTotal,sizeof(int),1,FilePoint);
	BeNumber=malloc(BlockTotal*sizeof(int));
	Beta=malloc(BlockTotal*sizeof(float));


	for(k=0;k<BlockTotal;k++)
	{	fread(&Beta[k],sizeof(float),1,FilePoint);
		fread(&BeNumber[k],sizeof(int),1,FilePoint);
	}


	fread(&BeginNumber,sizeof(int),1,FilePoint);
	fread(&SegmentTotal,sizeof(int),1,FilePoint);
	fread(&Rlarge,sizeof(float),1,FilePoint);
	fread(&LASTND,sizeof(int),1,FilePoint);
	fread(&LASTEL,sizeof(int),1,FilePoint);

	if((LASTND >LND)||(LASTEL>LEL))
	{printf("There isn’tenough memory!\n");exit(0);}


	Node=malloc(LASTND*sizeof(struct Nodetype));
	fread(Node,sizeof(struct Nodetype),LASTND,FilePoint);
	fread(&LASTEL,sizeof(int),1,FilePoint);

	Element=malloc(LASTEL*sizeof(struct elementtype));
	fread(Element,sizeof(struct elementtype),LASTEL,FilePoint);
	fclose(FilePoint);

	Address=malloc(LASTND*sizeof(int));

	KMatricStorageAddress(BeginNumber,LASTND,Address);
	nk=Address[LASTND -1]+1;
	KMatric=malloc(nk*sizeof(float));
	PMatric=malloc(LASTND*sizeof(float));

	for(k=0;k<BlockTotal;k++)
		Beta[k]=1/(Beta[k]*MU0);

	KPMatricSetup(Node,Address,Element,Beta,BeNumber,BlockTotal,KMatric,nk,PMatric,LASTND);

	/* Deal with boundary condition according to certain problem */

	{	int nbn,count=0;
		float DeltaTheta;
		struct boundaryNodetype *BoundaryNode;

		nbn=2*SegmentTotal+BeginNumber;			/* BeginNumber=1 */
		BoundaryNode=malloc(nbn*sizeof(struct boundaryNodetype));

		for(k=0;k<=SegmentTotal;k++)
		{count+=k;
		 BoundaryNode[k].bNode=count;
		 BoundaryNode[k].bptal=0.0;
		}

		for(k=1;k<=SegmentTotal;k++)
		{DeltaTheta=PI/(2*SegmentTotal);
		 BoundaryNode[k+SegmentTotal].bNode=(1+SegmentTotal)*SegmentTotal/2+k;
		 BoundaryNode[k+SegmentTotal].bptal=-Rlarge*sin(k*DeltaTheta);
		}

		BoundaryConditionDeal(BoundaryNode,nbn,KMatric,PMatric,Address,LASTND);
		free(BoundaryNode);
	}


	solveq(KMatric,nk,PMatric,Address,LASTND);
	free(BeNumber);
	free(Beta);
	free(KMatric);
	free(Address);


	/*Fine the coordinate ofisopotentialpoints*/

   {	    	int i,j,count,num=40,Number0,Number;
		int ei,ej,em;
		float Maxfi=0.0,Minfi=0.0,Deltafi,fi,p[200][2];

		if((FilePoint=fopen("emf2dout.dat","wb"))==NULL)
		{printf("cannotopen outputfile！ \n");exit(0);}

		fwrite(&num,sizeof(int),1,FilePoint);

		for(k=0;k<LASTND;k++)
		{if(PMatric[k]>Maxfi)Maxfi=PMatric[k];
		if(PMatric[k]<Maxfi)Minfi=PMatric[k];
		}

		Deltafi=(Maxfi-Minfi)/(num+1);

	for(i=1;i<=num;i++)
	{	fi=Minfi+Deltafi*i;
		count=0;
		Number0=0;
		Number=1;

		for(j=0;j<SegmentTotal;j++)
		{
			for(k=Number0;k<Number0+Number;k+=2)
 	 		{	ei=Element[k].i;
				ej=Element[k].j;
				em =Element[k].m;

				if((fi-PMatric[ei])*(fi-PMatric[ej])<=0.0)
				{p[count][0]=Node[ei].x+(fi-PMatric[ei])*(Node[ej].x-Node[ei].x)/(PMatric[ej]-PMatric[ei]);
				 p[count][1]=Node[ei].y+(fi-PMatric[ei])*(Node[ej].y-Node[ei].y)/(PMatric[ej]-PMatric[ei]);
				 count++;
				}

				if((fi-PMatric[ej])*(fi-PMatric[em])<=0.0)
				{p[count][0]=Node[ej].x+(fi-PMatric[ej])*(Node[em].x-Node[ej].x)/(PMatric[em]-PMatric[ej]);
				p[count][1]=Node[ej].y+(fi-PMatric[ej])*(Node[em].y-Node[ej].y)/(PMatric[em]-PMatric[ej]);
				count++;
				}

				if((fi-PMatric[em])*(fi-PMatric[ei])<=0.0)
				{p[count][0]=Node[em].x+(fi-PMatric[em])*(Node[ei].x-Node[em].x)/(PMatric[ei]-PMatric[em]);
				p[count][1]=Node[em].y+(fi-PMatric[em])*(Node[ei].y-Node[em].y)/(PMatric[ei]-PMatric[em]);
				count++;
				}
			}

			Number0+=Number;
			Number+=2;

		}

		mbbub(p,count);

		fwrite(&count,sizeof(int),1,FilePoint);

		fwrite(&fi,sizeof(float),1,FilePoint);

		for(k=0;k<count;k++)
		{	fwrite(&p[k][0],sizeof(float),1,FilePoint);
		 	fwrite(&p[k][1],sizeof(float),1,FilePoint);
		}
	}

   }


	fclose(FilePoint);

}



void KMatricStorageAddress(int BeginNumber,int LASTND,int *Address)

{	int k,count,width0,*width;
	width=(int*)malloc(LASTND*sizeof(int));
	width[0]=1;

	if(BeginNumber>=2)
		for(k=1;k<BeginNumber;k++)
			width[k]=2;

	width[BeginNumber]=BeginNumber+1;
	width0=BeginNumber+2;
	count=0;


	for(k=BeginNumber+1;k<LASTND;k++)
	{
	width[k]=width0;
	count++;
		if(count==width0-1)
		{width0++;
		count=0;
		}
	}

	Address[0]=0;

	for(k=1;k<LASTND -1;k++)
	Address[k]=Address[k-1]+width[k];

	free(width);
}



void KPMatricSetup(struct Nodetype *Node,int *Address,struct elementtype *Element,float *Beta,int *BeNumber,int nb,float *KMatric,int nk,float *PMatric,int np)


{	int i,k,k0=0,ads;float bi,bj,bm,ci,cj,cm,Delta,kr;

	for(i=0;i<nk;i++)
	KMatric[i]=0.0;

	for(i=0;i<np;i++)
	PMatric[i]=0.0;

	for(i=0;i<nb;i++)

	{for(k=k0;k<k0+BeNumber[i];k++)

	    {    bi=Node[Element[k].j].y-Node[Element[k].m].y;
		 bj=Node[Element[k].m].y-Node[Element[k].i].y;
		 bm=Node[Element[k].i].y-Node[Element[k].j].y;
		 ci=Node[Element[k].m].x-Node[Element[k].j].x;
		 cj=Node[Element[k].i].x-Node[Element[k].m].x;
		 cm=Node[Element[k].j].x-Node[Element[k].i].x;
		 Delta=(bi*cj-bj*ci)/2;kr=Beta[i]/4/Delta;

		 ads=Address[Element[k].i];
		 KMatric[ads]+=kr*(bi*bi+ci*ci);
		 ads=Address[Element[k].j];
		 KMatric[ads]+=kr*(bj*bj+cj*cj);
		 ads=Address[Element[k].m];
		 KMatric[ads]+=kr*(bm*bm +cm*cm);


		if(Element[k].i>Element[k].j)

		ads=Address[Element[k].i]-(Element[k].i-Element[k].j);
		else

		ads=Address[Element[k].j]-(Element[k].j-Element[k].i);

		KMatric[ads]+=kr*(bi*bj+ci*cj);

		if(Element[k].j>Element[k].m)
		ads=Address[Element[k].j]-(Element[k].j-Element[k].m);
		else
		ads=Address[Element[k].m]-(Element[k].m -Element[k].j);

		KMatric[ads]+=kr*(bj*bm +cj*cm);

		if(Element[k].m >Element[k].i)
		ads=Address[Element[k].m]-(Element[k].m -Element[k].i);
		else
		ads=Address[Element[k].i]-(Element[k].i-Element[k].m);

		KMatric[ads]+=kr*(bm*bi+cm*ci);

	    }

		k0+=BeNumber[i];
	}
}




void BoundaryConditionDeal(struct boundaryNodetype *BoundaryNode,int nbn,float *KMatric,float *PMatric,int *Address,int np)


{	int k,n,n1,nn1,n11,n12,n13;
	float uu,Kij;
	for(k=0;k<nbn;k++)

	{	n=BoundaryNode[k].bNode;
		uu=BoundaryNode[k].bptal;
		KMatric[Address[n]]=1.0;
		PMatric[n]=uu;

		if(n>0)
		{n11=n-(Address[n]-Address[n-1]-1);
			n12=n-1;
			for(n1=n11;n1<=n12;n1++)

			{
			 	PMatric[n1]=PMatric[n1]-uu*KMatric[Address[n]-(n-n1)];
				KMatric[Address[n]-(n-n1)]=0.0;
			}

		}

		n13=n+1;
		if(n13<=np-1)

		{for(n1=n13;n1<np;n1++)

			{	nn1=n1-n;
				Kij=0.0;

				if(nn1<=(Address[n1]-Address[n1-1]-1))
					Kij=KMatric[Address[n1]-nn1];

				PMatric[n1]=PMatric[n1]-uu*Kij;
				if(Kij!=0.0)KMatric[Address[n1]-nn1]=0.0;
			}
		}
	}
}






void mbbub(float (*p)[2],int n)

/*int n;float p[][2];*/

{	int m,k,i,j;
	float dx,dy;

	k=0;m =n-1;

	while(k<m)
	{	j=m -1;m =0;

		for(i=k;i<=j;i++)

			if(p[i][0]>p[i+1][0])
			{dx=p[i][0];p[i][0]=p[i+1][0];
			p[i+1][0]=dx;m =i;dy=p[i][1];
			p[i][1]=p[i+1][1];p[i+1][1]=dy;
			}

		j=k+1;
		k=0;

		for(i=m;i>=j;i--)
			if(p[i-1][0]>p[i][0])
			{dx=p[i][0];p[i][0]=p[i-1][0];
			p[i-1][0]=dx;k=i;dy=p[i][1];
			p[i][1]=p[i-1][1];p[i-1][1]=dy;
			}

	}
	return;
}