📄 shenliu.cpp
字号:
if((place_mat[j][2]==mp[k].material)&&(mp[k].B[0][0]<1.0e-20))
{
for(i=place_mat[j][0];i<=place_mat[j][1];i++)
{
pE[i].zhx=(pE[i].nd[0].x+pE[i].nd[1].x+pE[i].nd[2].x)/3;
pE[i].zhy=(pE[i].nd[0].y+pE[i].nd[1].y+pE[i].nd[2].y)/3;
for(int aa=0;aa<biank_z;aa++)
{
if((pE[i].zhx>=jb_ord[aa][0])&&(pE[i].zhx<=jb_ord[aa][2])&&(pE[i].zhy>=jb_ord[aa][1])&&(pE[i].zhy<=jb_ord[aa][3])) //有时间可以进行简单的修改，
{
if(((pE[i].zhx-guanjy_ord1)<=3*paijv)&&((pE[i].zhx-guanjy_ord1)>=1.0))
{
int shu1=0;
shu1=int((pE[i].zhx-guanjy_ord1-0.5)/paijv);
if((shu1==0)||(shu1==2))
{
pE[i].zhkx=gxshi_xishu[aa][0]*(paijv/2.0)+gxshi_xishu[aa][1];
pE[i].zhky=pE[i].zhkx;
}
else
{
pE[i].zhkx=gxshi_xishu[aa][0]*(fabs(pE[i].zhx-guanjy_ord1-0.5-(paijv/2.0)-shu1*paijv))+gxshi_xishu[aa][1];
pE[i].zhky=pE[i].zhkx;
}
for(l=0;l<3;l++)
{
for(m=0;m<3;m++)
{
pE[i].Aij[l][m]=( pE[i].b[l]*pE[i].b[m]*pE[i].zhkx +
pE[i].c[l]*pE[i].c[m]*pE[i].zhky)*pE[i].A;
}
}
}
if((pE[i].zhx-guanjy_ord1)>3*paijv)
{
int shu2=0;
shu2=int((pE[i].zhx-guanjy_ord1)/paijv);
pE[i].zhkx=gxshi_xishu[aa][0]*(fabs(pE[i].zhx-guanjy_ord1-shu2*paijv))+gxshi_xishu[aa][1];
pE[i].zhky=pE[i].zhkx;
for(l=0;l<3;l++)
{
for(m=0;m<3;m++)
{
pE[i].Aij[l][m]=( pE[i].b[l]*pE[i].b[m]*pE[i].zhkx +
pE[i].c[l]*pE[i].c[m]*pE[i].zhky)*pE[i].A;
}
}
}
if((pE[i].zhx-guanjy_ord1)<1.0)
{
pE[i].zhkx=gxshi_xishu[aa][0]*(fabs(pE[i].zhx-guanjy_ord1-1.0))+gxshi_xishu[aa][1];
pE[i].zhky=pE[i].zhkx;
for(l=0;l<3;l++)
{
for(m=0;m<3;m++)
{
pE[i].Aij[l][m]=( pE[i].b[l]*pE[i].b[m]*pE[i].zhkx +
pE[i].c[l]*pE[i].c[m]*pE[i].zhky)*pE[i].A;
}
}
}
}
}
}
}
}
}
printf("OK!\n");
//单元矩阵元素累加到总体矩阵相应的位置上////////////////////////////////////////////////////
printf("单元矩阵元素累加到总体矩阵相应的位置上...\n");
for(int idx=0;idx<element;idx++)
for(i=0;i<3;i++)
{
for(j=0;j<3;j++)
{pMatrix[ pE[idx].nd[i].number*iNode+pE[idx ].nd[j].number] += pE[idx].Aij[i][j];
}
pMf[ pE[idx].nd[i].number ]=0;
}
printf("OK!\n");
///////////边界条件的处理////////////
printf("边界条件的处理...\n");
double dBig=pow(10,20); //边界条件对角线扩大法处理所用的大数
for(k=0;k<num_border;k++)
{
pMatrix[bl[k].jd*iNode+bl[k].jd] += dBig;
pMf[bl[k].jd]=pMatrix[bl[k].jd*iNode+bl[k].jd]*bl[k].dx;
}
printf("OK!\n");
/////////////////////////////////////////////////////////////////////////////////
printf("调用全选主元高斯消去法函数解方程组...\n");
Gauss(pMatrix,pMf,iNode); //调用全选主元高斯消去法函数解方程组
// for(k=0;k<iNode;k++)
// printf("%d %lf\n",k,pMf[k]);
printf("OK!\n");
/////计算等势线坐标///////////////////////////////
printf("计算完成，开始等势线信息计算。\n");
// struct Dengqingkuang *ds=new Dengqingkuang[3*element];
ds=(Dengqingkuang*)malloc((3*element)*sizeof(Dengqingkuang));
int s=0;
int tiaoshu=(int(shuitou_zuigao-shuitou_zuidi))/dengshixian_cha-1;
double dengshixian[400];
for(k=0;k<3*element;k++) //初始化
{
for(l=0;l<tiaoshu;l++)
{
ds[k].deng[l].x1=0;
ds[k].deng[l].y1=0;
ds[k].deng[l].x2=0;
ds[k].deng[l].y2=0;
ds[k].deng[l].x3=0;
ds[k].deng[l].y3=0;
}
}
for (double ss=int(shuitou_zuigao)-dengshixian_cha;ss>shuitou_zuidi;ss -= dengshixian_cha)
{
dengshixian[s]=ss;
for (j=0;j<element;j++)
{
if(((pMf[pE[j].nd[0].number]>ss)&&(pMf[pE[j].nd[1].number]<ss))||((pMf[pE[j].nd[0].number]<ss)&&(pMf[pE[j].nd[1].number]>ss)))
{
ds[j].deng[s].x1=pE[j].nd[0].x-(pMf[pE[j].nd[0].number]-ss)/(pMf[pE[j].nd[0].number]-pMf[pE[j].nd[1].number])*(pE[j].nd[0].x-pE[j].nd[1].x);
ds[j].deng[s].y1=pE[j].nd[0].y-(pMf[pE[j].nd[0].number]-ss)/(pMf[pE[j].nd[0].number]-pMf[pE[j].nd[1].number])*(pE[j].nd[0].y-pE[j].nd[1].y);
}
if(((pMf[pE[j].nd[1].number]>ss)&&(pMf[pE[j].nd[2].number]<ss))||((pMf[pE[j].nd[1].number]<ss)&&(pMf[pE[j].nd[2].number]>ss)))
{
ds[j].deng[s].x2=pE[j].nd[1].x-(pMf[pE[j].nd[1].number]-ss)/(pMf[pE[j].nd[1].number]-pMf[pE[j].nd[2].number])*(pE[j].nd[1].x-pE[j].nd[2].x);
ds[j].deng[s].y2=pE[j].nd[1].y-(pMf[pE[j].nd[1].number]-ss)/(pMf[pE[j].nd[1].number]-pMf[pE[j].nd[2].number])*(pE[j].nd[1].y-pE[j].nd[2].y);
}
if(((pMf[pE[j].nd[2].number]>ss)&&(pMf[pE[j].nd[0].number]<ss))||((pMf[pE[j].nd[2].number]<ss)&&(pMf[pE[j].nd[0].number]>ss)))
{
ds[j].deng[s].x3=pE[j].nd[2].x-(pMf[pE[j].nd[2].number]-ss)/(pMf[pE[j].nd[2].number]-pMf[pE[j].nd[0].number])*(pE[j].nd[2].x-pE[j].nd[0].x);
ds[j].deng[s].y3=pE[j].nd[2].y-(pMf[pE[j].nd[2].number]-ss)/(pMf[pE[j].nd[2].number]-pMf[pE[j].nd[0].number])*(pE[j].nd[2].y-pE[j].nd[0].y);
}
}
s++;
}
printf("OK!\n");
/////////////////////////////////////////////////////
printf("写计算结果数据到文件...\n");
FILE *wfp; //文件操作
if((wfp=fopen("dat.txt","w+"))==NULL)
printf("Cann't open the file... ");
fprintf(wfp,"渗透系数为5E-4～E-5的情况。\n");
// fprintf(wfp,"%d\n",cishu);
fprintf(wfp,"等势线的坐标情况:\n");
for(j=0;j<s;j++)
{
fprintf(wfp,"%lf %d\n",dengshixian[j],-2);
/*/////最好按y值的大小进行排序//////////////////////////////
int array[SIZE]={99,46,4,80,30,2,8,68,5,37};
int i,j,t;
printf("原始数据是：\n");
for(i=0;i<=SIZE-1;i++)
printf("%4d",array[i]);
printf("\n");
for(j=0;j<=SIZE-1;j++)
for (i=0;i<=SIZE-2;i++)
if (array[i]>array[i+1])
{
t=array[i];
array[i]=array[i+1];
array[i+1]=t;
}
printf("按升序排列数据：\n");
for(i=0;i<=SIZE-1;i++)
printf("%4d",array[i]);
printf("\n");
*///////////////////////////////////////////////////
for(m=0;m<element;m++)
{
if((ds[m].deng[j].x1>1.0E-20)||(ds[m].deng[j].y1>1.0E-20))
fprintf(wfp,"%d %lf %lf\n",m,ds[m].deng[j].x1,ds[m].deng[j].y1);
if((ds[m].deng[j].x2>1.0E-20)||(ds[m].deng[j].y2>1.0E-20))
fprintf(wfp,"%d %lf %lf\n",m,ds[m].deng[j].x2,ds[m].deng[j].y2);
if((ds[m].deng[j].x3>1.0E-20)||(ds[m].deng[j].y3>1.0E-20))
fprintf(wfp,"%d %lf %lf\n",m,ds[m].deng[j].x3,ds[m].deng[j].y3);
}
}
fprintf(wfp,"计算得各结点上的水头值为:\n");
for(i=0;i<iNode;i++)
{
fprintf(wfp,"%d %f\n",i,pMf[i]);
}
printf("OK!\n");
fclose(wfp);
}
catch(...)
{
printf("Error occured...\n");
}
//释放总体矩阵和三角形单元数组占用内存
delete [] pMf; delete [] pMatrix; delete [] chf; delete [] bl;
delete [] mp; delete [] ds; // delete [] be;
free(pE); free(cnode);
printf("Please press Enter to exit...");
getchar();
}
//---------- 全选主元高斯消去法 ------------------------------
// a 体积为n*n的双精度实型二维数组，方程组系数矩阵，返回时将被破坏
// b 长度为n的双精度实型一维数组，方程组右端的常数向量，返回方程组的解向量
// n 整型变量，方程组的阶数
//--------------------------------------------------------------
int Gauss(double a[],double b[],int n)
{
int *js,l,k,i,j,is,p,q;
double d,t;
js=(int*)malloc(n*sizeof(int));
l=1;
for(k=0;k<=n-2;k++)
{
d=0.0;
for(i=k;i<=n-1;i++)
for(j=k;j<=n-1;j++)
{
t=fabs(a[i*n+j]);
if(t>d) { d=t; js[k]=j; is=i;}
}
if(d+1.0==1.0) l=0;
else
{
if(js[k]!=k)
for(i=0;i<=n-1;i++)
{
p=i*n+k; q=i*n+js[k];
t=a[p]; a[p]=a[q]; a[q]=t;
}
if(is!=k)
{
for(j=k;j<=n-1;j++)
{
p=k*n+j; q=is*n+j;
t=a[p]; a[p]=a[q]; a[q]=t;
}
t=b[k]; b[k]=b[is]; b[is]=t;
}
}
if(l==0)
{
free(js);
printf("Gauss funtion failed 1...\n");
return(0);
}
d=a[k*n+k];
for(j=k+1;j<=n-1;j++)
{
p=k*n+j; a[p]=a[p]/d;
}
b[k]=b[k]/d;
for(i=k+1;i<=n-1;i++)
{
for(j=k+1;j<=n-1;j++)
{
p=i*n+j;
a[p]=a[p]-a[i*n+k]*a[k*n+j];
}
b[i]=b[i]-a[i*n+k]*b[k];
}
}
d=a[(n-1)*n+n-1];
if(fabs(d)+1.0==1.0)
{
free(js);
printf("Gauss funtion failed 2...\n");
return(0);
}
b[n-1]=b[n-1]/d;
for(i=n-2;i>=0;i--)
{
t=0.0;
for(j=i+1;j<=n-1;j++)
t=t+a[i*n+j]*b[j];
b[i]=b[i]-t;
}
js[n-1]=n-1;
for(k=n-1;k>=0;k--)
if(js[k]!=k)
{
t=b[k]; b[k]=b[js[k]]; b[js[k]]=t;
}
free(js);
return(1);
}
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -