?? m.txt
字號:
# 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;
}?? 快捷鍵說明
復制代碼
Ctrl + C
搜索代碼
Ctrl + F
全屏模式
F11
切換主題
Ctrl + Shift + D
顯示快捷鍵
?
增大字號
Ctrl + =
減小字號
Ctrl + -