Пока готовил код для посторонних глаз, нашел косяк, исправил, теперь еще другая проблема. Решение расходится при размере ячейки меньше 1. Правда, этого можно избежать, если умножить коэффициенты в матрице на какое-то число, для шага 0,5, например, на 2, для 0,025 - на 600. (Этот момент там, в самом низу)
Вот код, немного большой, но я старался так сгруппировать, чтоб понятно было..
Код:
class vol
{
public:
unsigned char state; //byte 0 for w; 1 - e; 2 - n; 3 - s; байт, в котором хранится информация о том, какая сторона ячейки - граничная
coord a1; //координаты вершин
coord a2;
coord a3;
coord a4;
coord c; //координаты центра масс
float U;
float V;
float Vt;
float Ut;
float P;
int i;
int j;
vol(float x1, float y1,float x2, float y2,float x3, float y3,float x4, float y4,int i1,int j1)
{
a1.x = x1;
a1.y = y1;
a2.x = x2;
a2.y = y2;
a3.x = x3;
a3.y = y3;
a4.x = x4;
a4.y = y4;
i = i1;
j = j1;
U = 0;
V = 0;
Ut = 0;
Vt = 0;
P = 1;
state = 0;
}
void CompTemp(); //скорость без учета давления
void AddBToMatr(); //составление вектора правых частей для СЛАУ
void Correct();
void CompC()
{
float s1 = ((a1.x-a2.x)*(a4.y-a2.y)-(a4.x-a2.x)*(a1.y-a2.y))/2;
float s2 = -0.5*((a1.x-a3.x)*(a4.y-a3.y)-(a4.x-a3.x)*(a1.y-a3.y));
float v = 0.5*((a4.x-a1.x)*(a2.y-a3.y)-(a2.x-a3.x)*(a4.y-a1.y));
c.x = ((a1.x+a2.x+a4.x)/3*s1+(a1.x+a3.x+a4.x)/3*s2)/v; //сначала ищутся центры масс треугольников, потом
c.y = ((a1.y+a2.y+a4.y)/3*s1+(a1.y+a3.y+a4.y)/3*s2)/v; //по формуле - ц.м. четырехугольника
}
};
vector<vector<vol>> vec;
void vol::Correct()
{
float Pw = 0;
float Pe = 0;
float Pn = 0;
float Ps = 0;
P = x.items[(i)*(ny)+(j)]; //х - это вектор с давлением, который получаем в результате решения СЛАУ
if((bool)(state&1))
Pw = P;
else
Pw = x.items[(i-1)*(ny)+(j)];
if((bool)(state&2))
Pe = P;
else
Pe = x.items[(i+1)*(ny)+(j)];
if((bool)(state&4))
Pn = P;
else
Pn = x.items[(i)*(ny)+(1+j)];
if((bool)(state&8))
Ps = P;
else
Ps = x.items[(i)*(ny)+(j-1)];
float swy = (a2.x-a1.x); //компоненты площади
float ssy = -(a3.x-a1.x);
float sny = (a4.x-a2.x);
float sey = -(a4.x-a3.x);
float swx = -(a2.y-a1.y);
float ssx = (a3.y-a1.y);
float snx = -(a4.y-a2.y);
float sex = (a4.y-a3.y);
float v = 0.5*((a4.x-a1.x)*(a2.y-a3.y)-(a2.x-a3.x)*(a4.y-a1.y));
if(i == 1) //на входе и выходе градиент давления сохраняется
Pw = 2*P-Pe;
if(i == nx-2)
Pe = 2*P-Pw;
float pe = (Pe+P)/2;
float pw = (Pw+P)/2;
float ps = (Ps+P)/2;
float pn = (Pn+P)/2;
float t1 = (pe*sex+pn*snx+pw*swx+ps*ssx)*dt/(ro*v);
U = Ut-(pe*sex+pn*snx+pw*swx+ps*ssx)*dt/(ro*v);
V = Vt-(pe*sey+pn*sny+pw*swy+ps*ssy)*dt/(ro*v);
}
void vol::AddBToMatr()
{
float swy = (a2.x-a1.x);
float ssy = -(a3.x-a1.x);
float sny = (a4.x-a2.x);
float sey = -(a4.x-a3.x);
float swx = -(a2.y-a1.y);
float ssx = (a3.y-a1.y);
float snx = -(a4.y-a2.y);
float sex = (a4.y-a3.y);
float v = 0.5*((a4.x-a1.x)*(a2.y-a3.y)-(a2.x-a3.x)*(a4.y-a1.y));
float Ue = 0;
float Un = 0;
float Uw = 0;
float Us = 0;
float Ve = 0;
float Vn = 0;
float Vw = 0;
float Vs = 0;
Uw = vec[i-1][j].Ut;
Vw = vec[i-1][j].Vt;
Us = vec[i][j-1].Ut;
Vs = vec[i][j-1].Vt;
Ue = vec[i+1][j].Ut;
Ve = vec[i+1][j].Vt;
Un = vec[i][j+1].Ut;
Vn = vec[i][j+1].Vt;
float gw = (Ut+Uw)*swx+(Vt+Vw)*swy;
float ge = (Ut+Ue)*sex+(Vt+Ve)*sey;
float gn = (Ut+Un)*snx+(Vt+Vn)*sny;
float gs = (Ut+Us)*ssx+(Vt+Vs)*ssy;
if((bool)(state&128))
{
if((bool)(state&1))
gw = 0;
//на границах градиент давления равен нулю.
if((bool)(state&2))
ge = 0;
if((bool)(state&4))
gn = 0;
if((bool)(state&8))
gs = 0;
}
if(i == 0)
gw = UB*swx;
if(i == nx-1)
ge = UB*sex;
b.items[i*(ny)+j] = (ro/dt)*0.5*(ge+gn+gw+gs)/v;
}
void vol::CompTemp()
{
//иксовые и игриковые компоненты площади
float swy = (a2.x-a1.x);
float ssy = -(a3.x-a1.x);
float sny = (a4.x-a2.x);
float sey = -(a4.x-a3.x);
float swx = -(a2.y-a1.y);
float ssx = (a3.y-a1.y);
float snx = -(a4.y-a2.y);
float sex = (a4.y-a3.y);
float v = 0.5*((a4.x-a1.x)*(a2.y-a3.y)-(a2.x-a3.x)*(a4.y-a1.y));
float Ue = 0;
float Un = 0;
float Uw = 0;
float Us = 0;
float Ve = 0;
float Vn = 0;
float Vw = 0;
float Vs = 0;
float dxe = 0; //расстояние до соседней ячейки
float dxn = 0;
float dxw = 0;
float dxs = 0;
float dye = 0;
float dyn = 0;
float dyw = 0;
float dys = 0;
Uw = vec[i-1][j].U;
Vw = vec[i-1][j].V;
dxw = c.x-vec[i-1][j].c.x;
dyw = c.y-vec[i-1][j].c.y;
Us = vec[i][j-1].U;
Vs = vec[i][j-1].V;
dxs = c.x-vec[i][j-1].c.x;
dys = c.y-vec[i][j-1].c.y;
Ue = vec[i+1][j].U;
Ve = vec[i+1][j].V;
dxe = vec[i+1][j].c.x-c.x;
dye = vec[i+1][j].c.y-c.y;
Un = vec[i][j+1].U;
Vn = vec[i][j+1].V;
dxn = vec[i][j+1].c.x-c.x;
dyn = vec[i][j+1].c.y-c.y;
if(state&1)
{
dxw = 2*c.x-(a1.x+a2.x);
dyw = 2*c.y-(a1.y+a2.y);
}
if(state&2)
{
dxe = (a4.x+a3.x)-2*c.x;
dye = (a4.y+a3.y)-2*c.y;
}
if(state&4)
{
dxn = (a4.x+a2.x)-2*c.x;
dyn = (a4.y+a2.y)-2*c.y;
}
if(state&8)
{
dxs = 2*c.x-(a1.x+a3.x);
dys = 2*c.y-(a1.y+a3.y);
}
float mn = dxn*dxn+dyn*dyn;
float mw = dxw*dxw+dyw*dyw;
float ms = dxs*dxs+dys*dys;
float me = dxe*dxe+dye*dye;
float ve = 0.25*((V+Ve)*((V+Ve)*sey+(U+Ue)*sex));
float vn = 0.25*((V+Vn)*((V+Vn)*sny+(U+Un)*snx));
float vs = 0.25*((V+Vs)*((V+Vs)*ssy+(U+Us)*ssx));
float vw = 0.25*((V+Vw)*((V+Vw)*swy+(U+Uw)*swx));
float ue = 0.25*((U+Ue)*((V+Ve)*sey+(U+Ue)*sex));
float un = 0.25*((U+Un)*((V+Vn)*sny+(U+Un)*snx));
float us = 0.25*((U+Us)*((V+Vs)*ssy+(U+Us)*ssx));
float uw = 0.25*((U+Uw)*((V+Vw)*swy+(U+Uw)*swx));
if((bool)(state&128))
{
if((bool)(state&1)&(i>0&i<nx-1))
{
vw = 0;
uw = 0;
}
if((bool)(state&2))
{
ve = 0;
ue = 0;
}
if((bool)(state&4))
{
vn = 0;
un = 0;
}
if((bool)(state&8))
{
vs = 0;
us = 0;
}
}
float Ax;
float Ay;
Ax = (ue+un+us+uw)/v;
Ay = (ve+vn+vs+vw)/v;
float De = 0;
float Dw = 0;
float Ds = 0;
float Dn = 0;
if(sex!=0)
De = (sex*dxe/me);
if(swx!=0)
Dw = (swx*dxw/mw);
if(ssx!=0)
Ds = (ssx*dxs/ms);
if(snx!=0)
Dn = (snx*dxn/mn);
if(sey!=0)
De += (sey*dye/me);
if(swy!=0)
Dw += (swy*dyw/mw);
if(ssy!=0)
Ds += (ssy*dys/ms);
if(sny!=0)
Dn += (sny*dyn/mn);
if(state&1)
{
Uw = -U;
Vw = -V;
}
if(state&2)
{
Ue = -U;
Ve = -V;
}
if(state&4)
{
Un = -U;
Vn = -V;
}
if(state&8)
{
Us = -U;
Vs = -V;
}
float Dx = u*((U-Uw)*Dw+(Ue-U)*De+(U-Us)*Ds+(Un-U)*Dn)/v;
float Dy = u*((V-Vw)*Dw+(Ve-V)*De+(V-Vs)*Ds+(Vn-V)*Dn)/v;
Ut = U + (Dx-Ax)*dt;
Vt = V + (Dy-Ay)*dt;
float gw = 600*Dw; //Вот здесь! Если не умножать на 600, то метод расходится, для шага по сетке, меньше 1.
float gs = 600*Ds;
float gn = 600*Dn;
float ge = 600*De;
if((bool)(state&128))
{
if((bool)(state&1)&(i>0&i<nx-1))
gw = 0;
if((bool)(state&2))
ge = 0;
if((bool)(state&4))
gn = 0;
if((bool)(state&8))
gs = 0;
}
float a = (-gs); //это я в таком виде записываю матрицу
float b = (-gw);
float c = -(gn-gs-gw+ge);
float d = (ge);
float e = (gn);
MRow temp(a,b,c,d,e);
Matr[i*(ny)+j] = temp;
}
void cont() //Функция, которая должна все решать
{
for(int i = 1;i<nx-1;i++)
for(int j = 1;j<ny-1;j++)
vec[i][j].CompTemp();
for(int i = 1;i<nx-1;i++)
for(int j = 1;j<ny-1;j++)
vec[i][j].AddBToMatr();
BCG();
for(int i = 1;i<nx-1;i++)
for(int j = 1;j<ny-1;j++)
vec[i][j].Correct();
}
. Уж очень похожая картина по всем признакам.