📄 shenliu.cpp
字号:
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <iostream.h>
int Gauss(double a[],double b[],int n); //全选主元高斯消去法
struct ETNode{ //单元结点结构体
double x,y; //单元结点坐标
int number; //单元结点在总体区域划分中的编号
};
struct ElementTriangle{ //三角形单元结构体
ETNode nd[3]; //存储相对应的单元结点号
double zhx,zhy; //单元中心点坐标
double zhkx,zhky; //单元中心点的渗透系数
double a[3],b[3],c[3]; //基函数的系数
double A; //单元面积
double Aij[3][3]; //单元有限元特征式系数矩阵
};
struct Modulus_penetr{ //地层的渗透系数
int material; //地层编号
char name_mat[15]; //地层名称
double B[2][2]; //渗透系数
};
struct Char_infile{ //输入文件中的注释字符
char title[30];
};
struct Border{ //边界条件
int jd; //单元编号
double dx; //水头大小
};
//struct Border_ele{ //边界单元结构体
// int a,b,c,d;
//};
struct Dengord{ //等势线坐标
double x1,y1;
double x2,y2;
double x3,y3;
};
struct Dengqingkuang{ //等势线的水头
Dengord deng[400]; //等势线的条数400
};
//-------------- 全局变量 ---------------------------
int i,j,k,l,m, //循环指标
id; //单元的循环指标
// int iN[13]; //注意改变数组的大小(节点总数)
int Ele[6960]; //注意改变数组的大小(单元总数)
// double zhi[6]; //注意改变数组的大小(自由潜润面的节点总数)
// double tem[6]; //注意改变数组的大小(自由潜润面的节点总数)
// int sum_border_ele; //边界单元的各数
// int jinshui_point; //进水点节点编号
int num_border; //已知水头节点的总数
// int num_waterout; //已知潜润面的节点总数
// int outpoint[200]; //自由面节点
int sum_material; //地层渗透系数种类
int place_mat[200][3]; //不同地层分布位置
double D; //单元矩阵行列式值
int iNode,element; //节点数和单元数
int sum_fenqu; //不同渗透系数分区总数
double k1,k2; //渗透系数
int biank_z; //变渗透系数区域总数
double jb_ord[30][4]; //最多10个变渗透系数区域起点
double gxshi_xishu[30][2]; //对应区域内渗透系数表达式系数
double paijv; //灌浆孔排距
ElementTriangle* pE; //单元三角形结构体数组指针
Dengqingkuang* ds; //等势线数组指针
// double ai,bi,ci; //基函数的系数
// double water_error; //自由水面误差结束标准
struct ETNode *cnode; //存坐标
double shuitou_zuigao,shuitou_zuidi,guanjy_ord1;
double dengshixian_cha;
//-------------------- 主程序 -------------------------//
void main()
{
struct Char_infile *chf=new Char_infile[19]; //输入文件注释结构数组指针及分配内存
struct Border *bl=new Border[400]; //边界条件，结构数组指针及分配内存
FILE *cfptr;
if ((cfptr=fopen("shl_jb2.dat","r"))==NULL)
printf("无法打开文件.\n");
else
fscanf(cfptr,"%s",chf[0].title);
fscanf(cfptr,"%s",chf[1].title);
fscanf(cfptr,"%d",&iNode);
fscanf(cfptr,"%s",chf[2].title);
fscanf(cfptr,"%d",&element);
// fscanf(cfptr,"%s",chf[3].title);
// fscanf(cfptr,"%d",&sum_border_ele);
// struct Border_ele *be=new Border_ele[sum_border_ele];
// fscanf(cfptr,"%s",chf[4].title);
// for(i=0;i<sum_border_ele;i++)
// fscanf(cfptr,"%d %d %d %d",&be[i].a,&be[i].b,&be[i].c,&be[i].d);
// fscanf(cfptr,"%s",chf[5].title);
// fscanf(cfptr,"%d",&jinshui_point);
fscanf(cfptr,"%s",chf[3].title);
fscanf(cfptr,"%d",&num_border);
fscanf(cfptr,"%s",chf[4].title);
for(k=0;k<num_border;k++)
fscanf(cfptr,"%d %lf",&bl[k].jd,&bl[k].dx);
// fscanf(cfptr,"%s",chf[8].title);
// fscanf(cfptr,"%d",&num_waterout);
// fscanf(cfptr,"%s",chf[9].title);
// for(m=0;m<num_waterout;m++)
// fscanf(cfptr,"%d",&outpoint[m]);
//为总体矩阵,三角形单元数组,f函数向量等分配存储内存
pE=(ElementTriangle*)malloc(element*sizeof(ElementTriangle));
cnode=(ETNode*)malloc((iNode)*sizeof(ETNode));
double *pMatrix=new double[(iNode*iNode)*sizeof(double)]; //总体矩阵指针
double *pMf=new double[iNode*sizeof(double)]; //f向量指针
//初始化值为0,因为下面要累加总体矩阵
for(i=0;i<iNode*iNode;i++)
pMatrix[i]=0;
for(i=0;i<iNode;i++)
pMf[i]=0;
// fscanf(cfptr,"%s",chf[10].title);
// fscanf(cfptr,"%lf",&water_error);
fscanf(cfptr,"%s",chf[5].title);
for(i=0;i<element;i++)
fscanf(cfptr,"%d %d %d %d",&Ele[i],&pE[i].nd[0].number,&pE[i].nd[1].number,&pE[i].nd[2].number);
fscanf(cfptr,"%s",chf[6].title);
for(j=0;j<iNode;j++)
fscanf(cfptr,"%d %lf %lf",&cnode[j].number,&cnode[j].x,&cnode[j].y);
for(i=0;i<element;i++)
{
pE[i].nd[0].x=cnode[pE[i].nd[0].number].x;
pE[i].nd[0].y=cnode[pE[i].nd[0].number].y;
pE[i].nd[1].x=cnode[pE[i].nd[1].number].x;
pE[i].nd[1].y=cnode[pE[i].nd[1].number].y;
pE[i].nd[2].x=cnode[pE[i].nd[2].number].x;
pE[i].nd[2].y=cnode[pE[i].nd[2].number].y;
}
fscanf(cfptr,"%s",chf[7].title);
fscanf(cfptr,"%d",&sum_material);
fscanf(cfptr,"%s",chf[8].title);
struct Modulus_penetr *mp=new Modulus_penetr[sum_material]; //不同地层的渗透系数，结构数组指针及分配内存
for(id=0;id<sum_material;id++)
{
fscanf(cfptr,"%d %s",&mp[id].material,mp[id].name_mat);
fscanf(cfptr,"%lf %lf",&mp[id].B[0][0],&mp[id].B[0][1]);
fscanf(cfptr,"%lf %lf",&mp[id].B[1][0],&mp[id].B[1][1]);
}
fscanf(cfptr,"%s",chf[9].title);
fscanf(cfptr,"%d",&sum_fenqu);
fscanf(cfptr,"%s",chf[10].title);
for(l=0;l<sum_fenqu;l++)
fscanf(cfptr,"%d %d %d",&place_mat[l][0],&place_mat[l][1],&place_mat[l][2]);
fscanf(cfptr,"%s",chf[11].title);
fscanf(cfptr,"%lf",&shuitou_zuigao);
fscanf(cfptr,"%s",chf[12].title);
fscanf(cfptr,"%lf",&shuitou_zuidi);
fscanf(cfptr,"%s",chf[13].title);
fscanf(cfptr,"%lf",&dengshixian_cha);
fscanf(cfptr,"%s",chf[14].title);
fscanf(cfptr,"%d",&biank_z);
fscanf(cfptr,"%s",chf[15].title);
for(i=0;i<biank_z;i++)
fscanf(cfptr,"%lf %lf %lf %lf",&jb_ord[i][0],&jb_ord[i][1],&jb_ord[i][2],&jb_ord[i][3]);
fscanf(cfptr,"%s",chf[16].title);
for(i=0;i<biank_z;i++)
fscanf(cfptr,"%lf %lf",&gxshi_xishu[i][0],&gxshi_xishu[i][1]);
fscanf(cfptr,"%s",chf[17].title);
fscanf(cfptr,"%lf",&paijv);
fscanf(cfptr,"%s",chf[18].title);
fscanf(cfptr,"%lf",&guanjy_ord1);
fclose(cfptr);
/*
/// 输出读入的文件数据
printf("%s\n",chf[0].title);
printf("%s\n",chf[1].title);
printf("%d\n",iNode);
printf("%s\n",chf[2].title);
printf("%d\n",element);
// printf("%s\n",chf[3].title);
// printf("%d\n",sum_border_ele);
// printf("%s\n",chf[4].title);
// for(i=0;i<sum_border_ele;i++)
// printf("%d %d %d %d\n",be[i].a,be[i].b,be[i].c,be[i].d);
// printf("%s\n",chf[5].title);
// printf("%d\n",jinshui_point);
printf("%s\n",chf[3].title);
printf("%d\n",num_border);
printf("%s\n",chf[4].title);
for(k=0;k<num_border;k++)
printf("%d %lf\n",bl[k].jd,bl[k].dx);
// printf("%s\n",chf[8].title);
// printf("%d\n",num_waterout);
// printf("%s\n",chf[9].title);
// for(m=0;m<num_waterout;m++)
// printf("%d\n",outpoint[m]);
// printf("%s\n",chf[10].title);
// printf("%E\n",water_error);
printf("%s\n",chf[5].title);
for(i=0;i<element;i++)
printf("%d %d %d %d\n",Ele[i],pE[i].nd[0].number,pE[i].nd[1].number,pE[i].nd[2].number);
printf("%s\n",chf[6].title);
for(j=0;j<iNode;j++)
printf("%d %lf %lf\n",cnode[j].number,cnode[j].x,cnode[j].y);
printf("%s\n",chf[7].title);
printf("%d\n",sum_material);
printf("%s\n",chf[8].title);
for(id=0;id<sum_material;id++)
{
printf("%d %s\n",mp[id].material,mp[id].name_mat);
printf("%lf %lf\n",mp[id].B[0][0],mp[id].B[0][1]);
printf("%lf %lf\n",mp[id].B[1][0],mp[id].B[1][1]);
}
printf("%s\n",chf[9].title);
printf("%d\n",sum_fenqu);
printf("%s\n",chf[10].title);
for(l=0;l<sum_fenqu;l++)
printf("%d %d %d\n",place_mat[l][0],place_mat[l][1],place_mat[l][2]);
printf("%s\n",chf[11].title);
printf("%lf\n",shuitou_zuigao);
printf("%s\n",chf[12].title);
printf("%lf\n",shuitou_zuidi);
printf("%s\n",chf[13].title);
printf("%lf\n",dengshixian_cha);
printf("%s\n",chf[14].title);
printf("%d\n",biank_z);
printf("%s\n",chf[15].title);
for(i=0;i<biank_z;i++)
printf("%lf %lf %lf %lf\n",jb_ord[i][0],jb_ord[i][1],jb_ord[i][2],jb_ord[i][3]);
printf("%s\n",chf[16].title);
for(i=0;i<biank_z;i++)
printf("%lf %lf\n",gxshi_xishu[i][0],gxshi_xishu[i][1]);
printf("%s\n",chf[17].title);
printf("%lf\n",paijv);
printf("%s\n",chf[18].title);
printf("%lf\n",guanjy_ord1);
printf("节点信息.\n");
for(i=0;i<element;i++)
{
printf("%d %lf %lf\n",pE[i].nd[0].number,pE[i].nd[0].x,pE[i].nd[0].y);
printf("%d %lf %lf\n",pE[i].nd[1].number,pE[i].nd[1].x,pE[i].nd[1].y);
printf("%d %lf %lf\n\n",pE[i].nd[2].number,pE[i].nd[2].x,pE[i].nd[2].y);
}
*/
try{
printf("计算基函数系数值...\n");
for(id=0;id<element;id++)
{
for(i=0;i<3;i++)
{
if(i==0) j=1,k=2;
else if(i==1) j=2,k=0;
else if(i==2) j=0,k=1;
pE[id].A=((pE[id].nd[j].x-pE[id].nd[i].x)*(pE[id].nd[k].y-pE[id].nd[i].y)-
(pE[id].nd[j].y-pE[id].nd[i].y)*(pE[id].nd[k].x-pE[id].nd[i].x))/2.0;
D=2.0*pE[id].A;
pE[id].a[i]=( pE[id].nd[j].x*pE[id].nd[k].y- pE[id].nd[k].x*pE[id].nd[j].y )/D;
pE[id].b[i]=( pE[id].nd[j].y-pE[id].nd[k].y )/D;
pE[id].c[i]=( pE[id].nd[k].x-pE[id].nd[j].x )/D;
}
}
printf("OK!\n");
printf("计算单元渗透系数及有限元特征式系数矩阵...\n");
for(j=0;j<sum_fenqu;j++)
{
for(k=0;k<sum_material;k++)
{
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++)
{
k1=mp[k].B[0][0];
k2=mp[k].B[1][1];
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]*k1 +
pE[i].c[l]*pE[i].c[m]*k2)*pE[i].A;
}
}
}
}
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -