椭圆生成算法.docx
《椭圆生成算法.docx》由会员分享，可在线阅读，更多相关《椭圆生成算法.docx（96页珍藏版）》请在冰豆网上搜索。
椭圆生成算法
#include"EllipseGen.h"
structparameter0
{
doubleal,bl,cl,a,b,c,center_x,center_y,center_z,ymin,ymax,xmin,xmax,zmin,zmax;
doublealpha1,bata1,gama1,alpha2,bata2,gama2,alpha3,bata3,gama3;
doublegamamin,gamamax,miumin,miumax;
doubleA,B,C;
doublecoe[7];
};
structparameter
{
doubleal,bl,cl,a,b,c,center_x,center_y,center_z;
doubleA,B,C;
doublealpha1,bata1,gama1,alpha2,bata2,gama2,alpha3,bata3,gama3;
structparameter*next;
};
structbasicparameter
{
doubleal,bl,cl,a,b,c,center_x,center_y,center_z;
doublealpha1,bata1,gama1,alpha2,bata2,gama2,alpha3,bata3,gama3;
doubleA,B,C;
};
structcoordinate
{
doublex;
doubley;
};
charpercent='%';
doubleM=pow（2.0,42）;
doublelamda=pow（5.0,17）;
structparameter*head;
doubleb_aver=0,a_aver=0;
doublefamin=0,famax=0,fb=0;
doubleeamin=0,eamax=0,ebmin=0,ecmin=0;
doublecamin=0,camax=0,cbmin=0,cc=0;
doublegamamin=0,gamamax=0,miumin=0,miumax=0;
doublesitagama=0,miugama=0,sitamiu=0,miumiu=0;
doubleleftbar,rightbar,formbar,backbar,upbar,downbar,leftbar1,rightbar1,backbar1,formbar1,downbar1,upbar1;
intinitialx,finalx,initialy[2],finaly[2];
intinitialwz[2],finalwz[2],initialkz[2],finalkz[2];
double*total_xmin[2],*total_xmax[2],*total_ymin[2],*total_ymax[2];
doubleerror=pow（10.0,-8）,errort=pow（10.0,-8）;
doublef_xseed=SEED,f_yseed=SEED,f_zseed=SEED,f_aseed=SEED;
doublee_xseed=SEED,e_yseed=SEED,e_zseed=SEED,e_aseed=SEED,e_bseed=SEED,e_cseed=SEED;
doublec_xseed=SEED,c_yseed=SEED,c_zseed=SEED,c_aseed=SEED,c_bseed=SEED;
doubleig1seed=SEED,im1seed=SEED,ig2seed=SEED,im2seed=SEED;
doublemseed=SEED;
doubleseqx=SEED,seqy=SEED,seqz=SEED;
intEllipseGen:
:
uniell_generation（stringinput_file,stringout_file）
{
charrline[250];
intdistribution_sign1,distribution_sign2;
ifstreaminput_stream;
input_stream.open（input_file.c_str（）,ios:
:
in）;
if（!
input_stream）
{
hout<<"Cannotopenthefile:
"<<"\n";
return0;
};
istrstreamistr0（get_line（input_stream,rline））;
istr0>>leftbar>>rightbar>>backbar>>formbar
>>downbar>>upbar;//长方体区域的六条边
istrstreamistr2（get_line（input_stream,rline））;
istr2>>distribution_sign1;//倾角分布信息，取0为均匀分布，取1为正态分布
compact_elli（input_stream,distribution_sign1）;
istrstreamistr3（get_line（input_stream,rline））;
istr3>>distribution_sign2;//中心分布信息，取0为均匀分布，1为正态分布，2为指数分布
if（distribution_sign2==0）
uniform（input_stream,out_file）;
else
no_uniform（distribution_sign2,input_stream,out_file）;
return1;
};
intEllipseGen:
:
compact_elli（ifstream&input_stream,intdistribution_sign1）
{
inti,j,k,k1;
inti1,i2,j1;
intfm,em,cm,m;
inte_num=0,f_num=0,c_num=0;
intn,num=0;
intnorm1,norm2,norm4,norm5;
intf_e_c_division;//若要生成纤维取0，若要生成椭球取1，若要生成薄钱币颗粒取2
double*ixseed,*iyseed,*izseed;
doubledz,dy,dx;
charrline[250];
intshape_sign1,shape_sign2,shape_sign3;//shape1取1表示存在纤维颗粒，取0表示不存在
//同理，shape2,shape3取1代表椭球颗粒以及薄钱币颗粒存在，否则不存在
structparameter0**elliptic[3],*p0;
structparameter*p1,*p2;
/***********************确定有几种形态的颗粒，并读取相关的参数*************************/
istrstreamistr1（get_line（input_stream,rline））;
istr1>>shape_sign1>>shape_sign2>>shape_sign3;
/*输入纤维的aminamaxbmincmingamamingamamaxmiuminmiumax*/
//纤维长轴a服从（famin,famax）的均匀分布,中轴b取fb,短轴c与b同
//长轴a与z轴的夹角服从（fgamamin,fgamamax）的均匀分布,a在xoy面上投影与x轴的夹角服从（fmiumin,fmiumax）的均匀分布
if（shape_sign1==1）
{
istrstreamistr20（get_line（input_stream,rline））;
istr20>>famin>>famax>>fb;
}
/*输入椭球体的aminamaxbmincmingamamingamamaxmiuminmiumax*/
//椭球长轴a服从（eamin,eamax）的均匀分布,中轴b服从（ebmin,a）的均匀分布,短轴c服从（ecmin,b）的均匀分布
//长轴a与z轴的夹角服从（egamamin,egamamax）的均匀分布,a在xoy面上投影与x轴的夹角服从（emiumin,emiumax）的均匀分布
if（shape_sign2==1）
{
istrstreamistr21（get_line（input_stream,rline））;
istr21>>eamin>>eamax>>ebmin>>ecmin;
}
/****************************************/
/*输入薄钱币颗粒的aminamaxbmincmingamamingamamaxmiuminmiumax*/
//薄钱币颗粒长轴a服从（camin,camax）的均匀分布,中轴b服从（cbmin,a）,短轴c取cc
//长轴a与z轴的夹角服从（cgamamin,cgamamax）的均匀分布,a在xoy面上投影与x轴的夹角服从（cmiumin,cmiumax）的均匀分布
if（shape_sign3==1）
{
istrstreamistr22（get_line（input_stream,rline））;
istr22>>camin>>camax>>cbmin>>cc;
}
/*********************颗粒形态确定完毕**********************************************/
//颗粒方向参数的读取
if（distribution_sign1==0）
{
istrstreamistr30（get_line（input_stream,rline））;
istr30>>gamamin>>gamamax>>miumin>>miumax;
}
else
{
istrstreamistr31（get_line（input_stream,rline））;
istr31>>sitagama>>miugama>>sitamiu>>miumiu;
}
//纤维，椭球，薄钱币颗粒的个数比为fm:
em:
cm
istrstreamistr4（get_line（input_stream,rline））;
istr4>>fm>>em>>cm;
/*****************输入长方体区域参数**********************************/
m=fm+em+cm;
if（ebmin==0&&ecmin==0）
{
b_aver=eamin*em/（fm+em+cm）;
}
else
{
b_aver=（fb*fm/（fm+em+cm）+ebmin*em/（fm+em+cm）+cbmin*cm/（fm+em+cm））;//每一行可能容纳颗粒的最多个数由其决定
}
/*intsum=0;
if（shape_sign1==1）
{
sum=sum+3;
a_aver=a_aver+（famin+famax+fb）;
}
if（shape_sign2==1）
{
sum=sum+3;
a_aver=a_aver+（eamin+eamax+ebmin）;
}
if（shape_sign3==1）
{
sum=sum+3;
a_aver=a_aver+（camin+camax+cbmin）;
}
a_aver=a_aver/sum;*/
leftbar1=leftbar;
rightbar1=rightbar;
downbar1=downbar;
upbar1=upbar;
backbar1=backbar;
formbar1=formbar;
//申请空间
//elliptic[2]存储待生成层的颗粒，elliptic[0]，elliptic[1]存储此层前1层，前2层的颗粒
for（k=0;k<=2;k++）
{
elliptic[k]=（structparameter0**）malloc（N1*sizeof（structparameter0*））;
if（elliptic[k]==NULL）
{
printf（"Can'tobtaintheneededspace"）;
return0;
}
for（j=0;j{
elliptic[k][j]=（structparameter0*）malloc（N2*sizeof（structparameter0））;
if（elliptic[k][j]==NULL）
{
printf（"Can'tobtaintheneededspace"）;
return0;
}
}
}
for（k1=0;k1<=1;k1++）
{
total_xmin[k1]=（double*）malloc（N2*sizeof（double））;//存储待生成层的前1层，前2层每列x的最小值
total_xmax[k1]=（double*）malloc（N2*sizeof（double））;//存储待生成层的前1层，前2层每列x的最大值
total_ymin[k1]=（double*）malloc（N1*sizeof（double））;//存储待生成层的前1层，前2层每列y的最小值
total_ymax[k1]=（double*）malloc（N1*sizeof（double））;//存储待生成层的前1层，前2层每列y的最大值
}
//清零
for（k1=0;k1<=1;k1++）
for（i1=0;i1{
total_xmin[k1][i1]=0;
total_xmax[k1][i1]=0;
}
for（k1=0;k1<=1;k1++）
for（j1=0;j1{
total_ymin[k1][j1]=0;
total_ymax[k1][j1]=0;
}
for（k=0;k<=2;k++）
for（j=0;jfor（i=0;i{
elliptic[k][j][i].A=0;
elliptic[k][j][i].a=0;
elliptic[k][j][i].al=0;
elliptic[k][j][i].alpha1=0;
elliptic[k][j][i].alpha2=0;
elliptic[k][j][i].alpha3=0;
elliptic[k][j][i].B=0;
elliptic[k][j][i].b=0;
elliptic[k][j][i].bata1=0;
elliptic[k][j][i].bata2=0;
elliptic[k][j][i].bata3=0;
elliptic[k][j][i].bl=0;
elliptic[k][j][i].C=0;
elliptic[k][j][i].c=0;
elliptic[k][j][i].center_x=0;
elliptic[k][j][i].center_y=0;
elliptic[k][j][i].center_z=0;
elliptic[k][j][i].cl=0;
for（k1=0;k1<=6;k1++）
{
elliptic[k][j][i].coe[k1]=0;
}
elliptic[k][j][i].gama1=0;
elliptic[k][j][i].gama2=0;
elliptic[k][j][i].gama3=0;
elliptic[k][j][i].xmax=0;
elliptic[k][j][i].ymax=0;
elliptic[k][j][i].xmin=0;
elliptic[k][j][i].ymin=0;
elliptic[k][j][i].zmin=0;
elliptic[k][j][i].zmax=0;
}
n=0;
head=NULL;
//////******************************************椭球的定位******************************************************//////
for（k=0;;k++）
{
norm1=1;//颗粒的z最大值大于upbar1时为1，用于判断椭球生成是否完毕
for（j=0;;j++）
{
norm2=1;//颗粒y最大值大于formbar1时为1，用于判断此层椭球是否生成完毕
for（i=0;;i++）
{
norm4=1;//若椭球满足插入下一层的条件则取0
norm5=1;//若椭球满足插入下一行的条件则取0
//形成链表
{
n=n+1;
if（n==1）
{
p1=p2=（structparameter*）malloc（sizeof（structparameter））;
head=p1;
}
else
{
p1=（structparameter*）malloc（sizeof（structparameter））;
（*p2）.next=p1;
p2=p1;
}
}
p0=&elliptic[2][j][i];
doublem1;
m1=unifrnd（0,（double）（m）,&mseed）;
if（m1{
f_e_c_division=0;
f_num++;
}
elseif（m1>=fm&&m1<（em+fm））
{
f_e_c_division=1;
e_num++;
}
else
{
f_e_c_division=2;
c_num++;
}
//printf（"1j=%di=%dN1=%dN2=%d\n",j,i,N1,N2）;
parametergeneration1（p0,f_e_c_division,distribution_sign1）;
parametergeneration2（p0）;
if（f_e_c_division==0）
{
ixseed=&f_xseed;
iyseed=&f_yseed;
izseed=&f_zseed;
}
elseif（f_e_c_division==1）
{
ixseed=&e_xseed;
iyseed=&e_yseed;
izseed=&e_zseed;
}
else
{
ixseed=&c_xseed;
iyseed=&c_yseed;
izseed=&c_zseed;
}
///*********************************************x的第一次定位******************************************************///
if（i==0）
{//printf（"%lf%lf%lf%lf\n",cos（（*p0）.alpha1）,cos（（*p0）.alpha2）,（*p0）.alpha3,（*p0）.A）;
（*p0）.center_x=unifrnd（leftbar1,leftbar1+（a_aver+a_aver*h）,ixseed）;
（*p0）.xmin=（*p0）.center_x-sqrt（（*p0）.A）;
（*p0）.xmax=（*p0）.center_x+sqrt（（*p0）.A）;
}
else
{
for（i1=0;i1if（elliptic[2][j][i1].xmax>elliptic[2][j][i-1].xmin）
{
initialx=i1;
break;
}
finalx=i-1;//与当前椭球定位相关的当前行椭球，从第initialx个到第finalx个
（*p0）.center_x=elliptic[2][j][i-1].xmax+sqrt（（*p0）.A）;
（*p0）.xmin=（*p0）.center_x-sqrt（（*p0）.A）;
（*p0）.xmax=（*p0）.center_x+sqrt（（*p0）.A）;
}
///*********************************************x的第一次定位结束******************************************************///
////*********************************************y的第一次定位******************************************************///
if（j==0）
{
（*p0）.center_y=unifrnd（backbar1,backbar1+（a_aver+a_aver*h）,iyseed）;
（*p0）.ymin=（*p0）.center_y-sqrt（（*p0）.B）;
（*p0）.ymax=（*p0）.center_y+sqrt（（*p0）.B）;
}
else
{
//j为当前行，与当前颗粒定位相关的第j-
升级会员 