Расчет поля температур С++. Опции

[Krjuger]

В общем,я как всегда радую всех,новыми веселыми задачами.Сегодня на повеске дня следующая.
Рассчитать стационарное поле температур в пластине заданной формы. Предусмотреть возможность задания произвольных граничных условий (первого рода) и произвольной правой части (в коде программы).
ПОРЯДОК РЕШЕНИЯ ЗАДАЧИ:
1. Составить разностную схему “крест” для уравнения Пуассона в области заданной формы.
2. Проанализировать структуру получившейся системы линейных алгебраических уравнений и составить алгоритм расчета температурного поля.
3. Выписать расчетные формулы метода Зейделя и выбрать начальное приближение на основе задания граничного условия.
4. В качестве критерия окончания использовать условие . Предусмотреть подсчет числа итераций, потребовавшихся для достижения точности .
5. Решить задачу методом Зейделя.
6. Решить задачу методом релаксации. Параметр релаксации подобрать экспериментально.

Все фигуры составлены из квадратов с единичной стороной (т.о. вся фигура вписывается в квадрат 3х3). Начало координат выбирается в любом удобном для заданной фигуры месте.

Даааа,планы, как говорится грандиозные,не знаешь, как начинать разгребать.
Вот фото самой пластины.


К написанию кода я еще не приступил,нужна помощ в выборе реализации,а именно:
1)Как лучше реализовать хранение данных?одним массивом или как 7 массивов,которые с последствии будут склеиваться,чтоб выполнить пункт 1 задания,в общем,на данный момент я не совсем понимаю как лучше реализовать хранение данных,чтобы ими было удобно в последствии оперировать.

Буду рад любым идеям.

[Krjuger]

В общем я начал делать и пришел пока что к следующему.

#include<iostream>
#include<fstream>
#include<vector>
#include<math.h>
using namespace std;
const int N=2,m=3;
const double Eps=0.001;
vector<vector<double>> U,dU;

double Insert(double x,double y)
{
	return x*y;
}

double f(double x,double y)
{
	return x*(x*(y-x)+y);
}
double Cross(int i, int j)
{
	double temp=(1.0/4.0)*(U[i+1][j]+U[i-1][j]+U[i][j-1]+U[i][j+1]+(1.0/N*N)*f(double(j)/double(N),double(i)/double(N)));
	return temp;
}
double normi(vector<vector<double>> U)
{
	double S1=0;
	for(int i=1;i<U.size()-1;i++)
	{
		double S2=0;
		double n=min(U[i].size(),min(U[i-1].size(),U[i+1].size()));
		for(int j=1;j<n-1;j++)
		{
			S2+=abs(U[i][j]);
		}
		if (S2>S1) S1=S2;
	}
	return S1;
}
int main()
{
	freopen("info.txt","w",stdout);
	for(int i=0;i<=N;i++)
	{
		U.push_back(vector<double>(0));
		for(int j=N;j<=3*N;j++)
		{
			U[i].push_back(Insert(double(j)/double(N),double(i)/double(N)));
		}
	}
	for(int i=N+1;i<2*N;i++)
	{
		U.push_back(vector<double>(0));
		for(int j=0;j<=2*N;j++)
		{
			U[i].push_back(Insert(double(j)/double(N),double(i)/double(N)));
        }
	}
	for(int i=2*N;i<=3*N;i++)
	{
		U.push_back(vector<double>(0));
		for(int j=0;j<=3*N;j++)
		{
			U[i].push_back(Insert(double(j)/double(N),double(i)/double(N)));
		}
	}
	cout<<'\n';
	for(int i=0;i<U.size();i++)
	{
		for(int j=0;j<=(U[i].size()-1);j++)
		{
			cout<<U[i][j]<<"\t ";
		}
		cout<<'\n';
	}
	fclose(stdout);
	
	dU=U;
	double norm=normi(dU);
	int k=0;
	for(;norm>Eps;k++)
	{
		
		for(int i=1;i<U.size()-1;i++)
		{
			double n=min(U[i].size(),min(U[i-1].size(),U[i+1].size()));
			for(int j=1;j<n-1;j++)
			{
				dU[i][j]=U[i][j];
				U[i][j]=Cross(i,j);
				dU[i][j]=abs(dU[i][j]-U[i][j]);
			}
		}
		norm=normi(dU);
		freopen("info.txt","a",stdout);
	cout<<'\n';
	for(int i=0;i<U.size();i++)
	{
		for(int j=0;j<U[i].size();j++)
		{
			cout<<U[i][j]<<"\t ";
		}
		cout<<'\n';
	}
	cout<<"\n"<<"kol-vo iteraciy="<<k; 
	fclose(stdout);
	}
	return 0;
}
}

Так как не важно, как можно повернуть фигуру и с какой угол(сторону) считать точкой отсчета,я развернул фигуру на 90 градусов против часовой оси и за начало взял верхний левый угол.
Но тут возникли некоторые проблемы, если N изменить на 2 или более,то начинается пересчет точек,которые считаться не должны,т.к не хватает 1 данного.

                                  0	 0	 0	 0	 0	 
                                  0.5	 0.75	 1	 1.25	 1.5	 
                                  1	 1.5	 2	 2.5	 3	 
0          0.75	 1.5	 2.25	 3	 
0          1                   	 2	 3	 4	 5	 6	 
0          1.25	 2.5	 3.75	 5	 6.25	 7.5	 
0          1.5 	 3	 4.5	 6	 7.5	 9	 

Исходная матрица.
                                       0	                0	0	 0	 0	 
                                       0.5	                 0.8125	1.01563	 0.878906	 1.5	 
                                       1	                 1.29688	 1.82813	 2.08301	 3	 
0	1.19922	      2.31934	 3.16309	 3	 
0	1.45605	      3.06885	 4.52673	 4	 5              6	 
0	1.80151	      3.90509	 5.98296	 7.30824	 8.38956	 7.5	 
0	1.5	      3	                 4.5	 6	 7.5  	 9
Матрица после первой итерации.

Как видно элементы 0.75(нижний левый) и 2.5(верхний правый) не должны пересчитываться так,как не хватает по 1 числу.Да и в общем то я не совсем уверен,что представление данных в памати компьютера верное и в уравнение Пуассона подставляются именно необходимые элементы.Руками я подставлял и считал для матрици с N=1,там вроде верно,а дальше я не уверен,что будет правильно при других N.

Сообщение отредактировано: Krjuger - 2.05.2011 15:02