二维有限体积法计算热传导及源码

cout<<endl<<"不?能¨¹打䨰开a这a个?文?件t"<<str1<<endl; exit(1); } cout<<"load the data from "<<str1<<endl; fin>>L>>m>>n>>T; cout<<endl<<"The data loaded is:"<<endl;
nn=maxtrix[i][kk]/maxtrix[kk][kk]; b[i]=b[i]-nn*b[kk]; for(j=kk+1;j<mm;j++) {
maxtrix[i][j]=maxtrix[i][j]-nn*maxtrix[kk][j]; } i++; } kk++; } x[mm-1]=b[mm-1]/maxtrix[mm-1][mm-1]; for(kk=mm-2;kk>=0;kk--) { d=0.0; for(j=kk+1;j<mm;j++) {
if(iห้องสมุดไป่ตู้n)
{
if(i==0)
{
maxtrix[i][i+1]=-an[i];
maxtrix[i][i+n]=-ae[i];
}
else if(i==n-1)
{
maxtrix[i][i-1]=-as[i];
maxtrix[i][i+n]=-ae[i];
}
else
{
maxtrix[i][i+1]=-an[i];
maxtrix[i][i-1]=-as[i];
maxtrix[i][i+n]=-ae[i]; } } else if(i>=mm-n) { if(i==mm-1) {
maxtrix[i][i-1]=-as[i]; maxtrix[i][i-n]=-aw[i]; }
else if(i==mm-n) {
maxtrix[i]=new double [mm]; for(i=0;i<mm;i++)
b[i]=su[i]; for(i=0;i<mm;i++)
for(j=0;j<mm;j++) maxtrix[i][j]=0.0;
for(i=0;i<mm;i++) {
for(j=0;j<mm;j++)
{
maxtrix[i][i]=ap[i];
//#include "stdafx.h" #include <iostream> #include <vector> #include <math.h>
#include <cmath> #include <fstream> using namespace std; #define q 500 #define k 1000 void main() { //input the value you want:
sp[i]=0.0; ap[i]=aw[i]+ae[i]+as[i]+an[i]-sp[i]; } } //....................../// double *b=new double[mm];
double **maxtrix; maxtrix=new double *[mm]; for(i=0;i<mm;i++)
{ if(i<n-1) //第̨²一°?列¢D结¨¢点Ì? { aw[i]=0.0; ae[i]=k*Ax/dx; an[i]=k*Ay/dy; as[i]=k*Ay/dy; su[i]=q*Ax*k; sp[i]=0.0; ap[i]=aw[i]+ae[i]+as[i]+an[i]-sp[i]; } else if(i==n-1) //左Á¨®上¦?角?顶£¤点Ì? { as[i]=k*Ay/dy; ae[i]=k*Ax/dx; aw[i]=0.0; an[i]=0.0; sp[i]=-2*as[i]; su[i]=q*Ax*k+2*as[i]*T; ap[i]=aw[i]+ae[i]+as[i]+an[i]-sp[i]; } else if(i==mm-n) //右®¨°下?角?顶£¤点Ì? { ae[i]=0.0; as[i]=0.0; aw[i]=k*Ax/dx; an[i]=k*Ay/dy; su[i]=0; sp[i]=0; ap[i]=aw[i]+ae[i]+as[i]+an[i]-sp[i]; }
maxtrix[i][i+1]=-an[i];
maxtrix[i][i-n]=-aw[i]; } else {
maxtrix[i][i+1]=-an[i]; maxtrix[i][i-n]=-aw[i]; maxtrix[i][i-1]=-as[i]; } } else if(i%n==0) { maxtrix[i][i+1]=-an[i]; maxtrix[i][i+n]=-ae[i]; maxtrix[i][i-n]=-aw[i]; } else if((i+1)%n==0) { maxtrix[i][i-1]=-as[i]; maxtrix[i][i+n]=-ae[i]; maxtrix[i][i-n]=-aw[i]; } else { maxtrix[i][i+1]=-an[i]; maxtrix[i][i+n]=-ae[i]; maxtrix[i][i-n]=-aw[i]; maxtrix[i][i-1]=-as[i]; } } } cout<<endl; //............................for test.............// cout<<"The right const array is: "<<endl; for(i=0;i<mm;i++) cout<<b[i]<<" "; cout<<endl<<endl; cout<<"The maxtrix is :"<<endl; for(i=0;i<mm;i++) { for(j=0;j<mm;j++) cout<<maxtrix[i][j]<<" "; cout<<endl;
else if(i==mm-1) //右®¨°上¦?角?顶£¤点Ì? {
ae[i]=0.0; aw[i]=k*Ax/dx; as[i]=k*Ay/dy; an[i]=0.0; su[i]=2*as[i]*T; sp[i]=-2*as[i]; ap[i]=aw[i]+ae[i]+as[i]+an[i]-sp[i]; } else if(i%n==0 ) //最Á?底Ì¡Á下?一°?行D {
} //............................for test.............// //solve the maxtrix double *x,nn;
x=new double[mm]; kk=0; while(kk<mm-1) { i=kk+1; while(i<mm) {
double L,dx,dy,T,Ax,Ay,d; int m,n,i,j,kk,mm; //char str1[20]; char *str1="conf.dat"; //cout<<"请?输º?入¨?文?件t名?：êo"<<endl; cin>>str1; //string str1("data."); ifstream fin(str1); if(!fin) {
as[i]=0.0; ae[i]=k*Ax/dx; aw[i]=k*Ax/dx; an[i]=k*Ay/dy; sp[i]=0.0; su[i]=0.0; ap[i]=aw[i]+ae[i]+as[i]+an[i]-sp[i]; } else if((i+1)%n==0) //最Á?顶£¤上¦?一°?行D { an[i]=0.0; ae[i]=k*Ax/dx; aw[i]=k*Ax/dx; as[i]=k*Ay/dy; su[i]=2*as[i]*T; sp[i]=-2*as[i]; ap[i]=aw[i]+ae[i]+as[i]+an[i]-sp[i]; } else if(i>mm-n) //最Á?后¨®一°?列¢D { ae[i]=0; aw[i]=k*Ax/dx; an[i]=k*Ay/dy; as[i]=an[i]; su[i]=0.0; sp[i]=0.0; ap[i]=aw[i]+ae[i]+as[i]+an[i]-sp[i]; } else //内¨²点Ì? { ae[i]=k*Ax/dx; aw[i]=k*Ax/dx; an[i]=k*Ay/dy; as[i]=an[i]; su[i]=0.0;
d=d+maxtrix[kk][j]*x[j]; } x[kk]=(b[kk]-d)/maxtrix[kk][kk]; } //..............for test.............// cout<<endl<<"The result is: "<<endl; for(i=0;i<mm;i++) cout<<x[i]<<" "; cout<<endl; delete[] b,x; for(i=0;i<mm;i++) delete[] maxtrix[i]; delete[] maxtrix; }
cout<<"the thickness L="<<L<<"and T="<<T<<endl; cout<<"x axis has m="<<m<<" nodes, y axis has n="<<n<<" nodes"<<endl; //creat the maxtrix mm=m*n; dx=0.3/m; dy=0.4/n; vector<double> ae(mm),ap(mm),aw(mm),an(mm),as(mm),su(mm),sp(mm); Ax=L*dy; Ay=L*dx; //结¨¢点Ì?0,左Á¨®下?角?顶£¤点Ì? aw[0]=0.0; ae[0]=k*Ax/dx;an[0]=k*Ay/dy;as[0]=0.0; sp[0]=0.0;su[0]=q*Ax*k; ap[0]=aw[0]+ae[0]+an[0]+as[0]-sp[0]; for(i=1;i<mm;i++)