?? 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 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 ETNoded{ //迭代增加的單元結點結構體
// double x,y; //單元結點坐標
// int number; //單元結點在總體區域劃分中的編號
//};
struct Border_ele{ //邊界單元結構體
int a,b,c,d;
};
//-------------- 全局變量 ---------------------------
int i,j,k,l,m, //循環指標
id; //單元的循環指標
// int iN[13]; //注意改變數組的大小(節點總數)
int Ele[19]; //注意改變數組的大小(單元總數)
double zhi[16]; //注意改變數組的大小(節點總數)
double tem[16]; //注意改變數組的大小(節點總數)
// int tem1[4]; //注意改變數組的大小(可能的出水節點總數)
// int sum_border_ele; //邊界單元的各數
// int jinshui_point; //進水點節點編號
int num_border; //已知水頭節點的總數
int num_waterout; //可能出水點的節點總數
int outpoint[200]; //可能出水節點
int sum_material; //地層滲透系數種類
int place_mat[40][3]; //不同地層分布位置
double D; //單元矩陣行列式值
int iNode,element; //節點數和單元數
// ElementTriangle* pE; //單元三角形結構體數組指針
double ai,bi,ci; //基函數的系數
double water_error; //自由水面誤差結束標準
//-------------------- 主程序 -------------------------//
void main()
{
struct Char_infile *chf=new Char_infile[13]; //輸入文件注釋結構數組指針及分配內存
struct Border *bl=new Border[400]; //邊界條件，結構數組指針及分配內存
FILE *cfptr;
if ((cfptr=fopen("shu.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);
//為總體矩陣,三角形單元數組,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向量指針
struct ElementTriangle *pE=new ElementTriangle[element]; //單元
struct ETNode *cnode=new ETNode[iNode]; //存坐標
////"50"為認為規定迭代次數，要求計算次數較多的，根據情況自行改變此值//
//初始化值為0,因為下面要累加總體矩陣
for(i=0;i<iNode*iNode;i++)
pMatrix[i]=0;
for(i=0;i<iNode;i++)
pMf[i]=0;
fscanf(cfptr,"%s",chf[3].title);
fscanf(cfptr,"%lf",&water_error);
fscanf(cfptr,"%s",chf[4].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[5].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[6].title);
fscanf(cfptr,"%d",&num_border);
fscanf(cfptr,"%s",chf[7].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]);
fscanf(cfptr,"%s",chf[10].title);
fscanf(cfptr,"%d",&sum_material);
fscanf(cfptr,"%s",chf[11].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[12].title);
for(l=0;l<sum_material;l++)
fscanf(cfptr,"%d %d %d",&place_mat[l][0],&place_mat[l][1],&place_mat[l][2]);
/// 輸出讀入的文件數據
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("%E\n",water_error);
printf("%s\n",chf[4].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[5].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[6].title);
printf("%d\n",num_border);
printf("%s\n",chf[7].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("%d\n",sum_material);
printf("%s\n",chf[11].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[12].title);
for(l=0;l<sum_material;l++)
printf("%d %d %d\n",place_mat[l][0],place_mat[l][1],place_mat[l][2]);
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_material;j++)
{
for(i=place_mat[j][0];i<=place_mat[j][1];i++)
{
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]*mp[j].B[0][0] +
pE[i].c[l]*pE[i].c[m]*mp[j].B[1][1])*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,25); //邊界條件對角線擴大法處理所用的大數
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(i=0;i<iNode;i++)
?? 快捷鍵說明
復制代碼
Ctrl + C
搜索代碼
Ctrl + F
全屏模式
F11
切換主題
Ctrl + Shift + D
顯示快捷鍵
?
增大字號
Ctrl + =
減小字號
Ctrl + -