Решается задача теплопроводности для треугольной пластины: (0;0) (8;0) (4;4). На косых границах температура задана 100, на горизонтальной 50.
Результаты выводятся в файл, предназначенный для просмотра в MeshViewer.
Нужно сделать в точке на косой границе(например в точке (6;2) ) граничное условие третьего рода. Как это реализовать?)
Код:
// Implicit_Sheme_XYT.cpp: определяет точку входа для консольного приложения.
//
#include <math.h>
#include <fstream>
#include <iostream>
using namespace std;
double xc=4,xmax=8,ymax=4;
int K;
int M;
double t1=50;
double t2=100;
double t0=0;
int k1(int m)
{
return m;
}
int k2(int m)
{
return K-m;
}
int m2(int k)
{
return (k<=K/2?k:M-k+K/2);
}
int main()
{
int k,n,m;
double A=0.1; // du/dt=A*(d2u/dx2+d2u/dy2) A = lambda/c*ro
double u1=50;
double u2=97.72;
double U0=0;
double h=0.1;
K=(int)(xmax/h); //Количество откезков по x
M=(int)(ymax/h); //Количество откезков по y
double tau=0.1;
int N=(int)(250/tau); //Количество откезков по t
double R=tau*A/(h*h)/2;
double ***u = new double**[N+1];
for(n=0;n<=N;n++)
{
u[n] = new double*[K+1];
for(k=0;k<=K;k++)
u[n][k] = new double[m2(k)+1];
}
for (k=0;k<=K;k++)
{
u[0][k][0]=u1;
u[0][k][m2(k)]=u2;
}
/*for (k=0;k<=K;k++)
{
if (m2(k)==60)
{
u[0][k][m2(k)]=100*sqrt(2*h*h)+u[0][k-1][m2(k)-1];
}
else
{
u[0][k][m2(k)]=u2;
}
}*/
for (k=1;k<=K-1;k++)
for (m=1;m<=m2(k)-1;m++)
u[0][k][m]=U0;
double **u_ = new double*[K+1];
for(k=0;k<=K;k++)
u_[k] = new double[m2(k)+1];
//метод прогонки
double *L1 = new double[K+1];
double **M1 = new double*[K+1];
for(k=0;k<=K;k++)
M1[k] = new double[m2(k)+1];
double *L2 = new double[M+1];
double **M2 = new double*[K+1];
for(k=0;k<=K;k++)
M2[k] = new double[m2(k)+1];
L1[0]=0;
for (k=0;k<=K/2;k++)
M1[k][m2(k)]=u2;
for (k=0;k<=K-1;k++)
L1[k+1]=R/(-R*L1[k]+1+2*R);
L2[0]=0;
for (k=0;k<=K;k++)
M2[k][0]=u1;
for (m=0;m<=M-1;m++)
L2[m+1]=R/(-R*L2[m]+1+2*R);
for (n=0;n<=N-1;n++)
{
//Прямая
for (k=0;k<=K-1;k++)
for (m=0;m<=m2(k);m++)
M1[k+1][m]=(u[n][k][m]+R*M1[k][m])/(-R*L1[k]+1+2*R);
//Обратная
for (k=K/2;k<=K;k++)
u_[k][m2(k)]=u2;
for (k=K-1;k>=0;k--)
for (m=0;m<=m2(k)-1;m++)
u_[k][m]=L1[k+1]*u_[k+1][m]+M1[k+1][m];
for (k=0;k<=K;k++)
for (m=0;m<=m2(k)-1;m++)
M2[k][m+1]=(u_[k][m]+R*M2[k][m])/(-R*L2[m]+1+2*R);
//О
for (k=0;k<=K;k++)
u[n+1][k][m2(k)]=u2;
for (k=0;k<=K;k++)
for (m=m2(k)-1;m>=0;m--)
u[n+1][k][m]=L2[m+1]*u[n+1][k][m+1]+M2[k][m+1];
}
/*
if (m2(k)=60)
{
u[0][k][m2(k)]=100*sqrt(2*h*h)+u[0][k-1][m2(k)-1];
}
else
{
u[0][k][m2(k)]=u2;
}
if (m2(k)=60)
{
M1[k][m2(k)]=100*sqrt(2*h*h)+M1[k-1][m2(k)-1];
}
else
{
M1[k][m2(k)]=u2;
}
*/
std::ofstream fout("triangle.mv");
fout<<(K+1)*(M+1)<<" 3 7 U(0) U(1) U(10) U(25) U(250) \n";
for (k=0;k<=K;k++)
{
for (m=0;m<=M;m++)
if( m < m2(k))
fout <<k*(M+1)+m+1<<" "<<h*k<<" "<<h*m<<" 0 "<<u[0][k][m]<<' '<<u[1][k][m]<<' '<<u[N/250][k][m]<<' '<<u[N/10][k][m]<<' '<<u[N][k][m]<<' '<<h*k<<' '<<h*m<<"\n";
else fout <<k*(M+1)+m+1<<" "<<h*k<<" "<<h*m<<" 0 "<<u2+2.28<<' '<<u2+2.28<<' '<<u2+2.28<<' '<<u2+2.28<<' '<<u2+2.28<<" "<<h*k<<" "<<h*m<<"\n";
}
int node=0;
for (k=0;k<=K-1;k++)
for (m=0;m<=M-1;m++)
if(k>=k1(m)&&k<=k2(m)-1) node++;
fout <<node<<" 4 0\n";
int node1=0;
for (k=0;k<=K-1;k++)
for (m=0;m<=M-1;m++)
if(k>=k1(m)&&k<=k2(m)-1)
if((m==m2(k)-1)&&(k==k2(m)-1))
fout <<(++node1)<<" "<<k*(M+1)+m+1 <<" "<<(k+1)*(M+1)+m+1<<" "<<k*(M+1)+(m+1)+1<<" "<<k*(M+1)+(m+1)+1<<"\n";
else if((m==m2(k))&&(k==k1(m)))
fout <<(++node1)<<" "<<k*(M+1)+m+1<<" "<<(k+1)*(M+1)+m+1<<" "<<(k+1)*(M+1)+(m+1)+1<<" "<<(k+1)*(M+1)+(m+1)+1<<"\n";
else fout <<(++node1)<<" "<<k*(M+1)+m+1 <<" "<<(k+1)*(M+1)+m+1<<" "<<(k+1)*(M+1)+(m+1)+1<<" "<<k*(M+1)+(m+1)+1<<"\n";
for(n=0;n<N+1;n++)
{
for(k=0;k<K+1;k++)
delete[]u[n][k];
delete[]u[n];
}
delete[]u;
for(k=0;k<K+1;k++)
delete[]u_[k];
delete[]u_;
delete[]L1;
delete[]L2;
for(k=0;k<K+1;k++)
delete[]M2[k];
delete[]M2;
return 0;
}
20.06.2011, 11:45 
