# Разработка программы решения гармонической задачи методом конечных элементов (Лабораторная работа № 3), страница 2

{

slau.n = nx * ny * 2;

if(!(slau.ig = new int[slau.n + 1]))return 4;

generation_ig();

if(!(slau.jg = new int[slau.ig[slau.n]]))return 5;

generation_jg();

if(!(slau.ggl = new double[slau.ig[slau.n]]))return 6;

if(!(slau.ggu = new double[slau.ig[slau.n]]))return 7;

if(!(slau.di  = new double[slau.n]))return 8;

if(!(slau.pr  = new double[slau.n]))return 9;

return 0;

}

void MKE::generation_ig()

{

int i, j;

slau.ig[0] = 0;

slau.ig[1] = 0;

slau.ig[2] = 1;

for(i = 3; i < nx * 2; i += 2)

{

slau.ig[i]    = slau.ig[i-1] + 2;

slau.ig[i+1] = slau.ig[i]   + 3;

}

while(i <= slau.n)

{

slau.ig[i]     = slau.ig[i - 1] + 4;

slau.ig[i + 1] = slau.ig[i]     + 5;

for(j = 1, i += 2; j < nx - 1; j++, i += 2)

{

slau.ig[i]     = slau.ig[i - 1] + 8;

slau.ig[i + 1] = slau.ig[i]     + 9;

}

slau.ig[i]     = slau.ig[i - 1] + 6;

slau.ig[i + 1] = slau.ig[i]     + 7;

i += 2;

}

}

void MKE::generation_jg()

{

int i, j, k, t, nx2 = nx * 2;

slau.jg[0] = 0;

for(i = 2, j = 1; i < nx2; i += 2, j += 5)

{

slau.jg[j + 2] = slau.jg[j    ] = i - 2;

slau.jg[j + 3] = slau.jg[j + 1] = i - 1;

slau.jg[j + 4] = i;

}

while(i < slau.n)

{

for(k = 0; k < 4; k++, j++)

slau.jg[j + 4] = slau.jg[j] = i - nx2 + k;

slau.jg[j + 4] = i;

for(t = 2, i += 2, j += 5; t < nx; t++, i += 2, j += 10)

{

for(k = -2; k < 4; k++, j++)

slau.jg[j + 8] = slau.jg[j] = i - nx2 + k;

slau.jg[j + 8] = slau.jg[j] = i - 2;

j++;

slau.jg[j + 8] = slau.jg[j] = i - 1;

slau.jg[j + 9] = i;

}

for(k = -2; k < 2; k++, j++)

slau.jg[j + 6] = slau.jg[j] = i - nx2 + k;

slau.jg[j + 6] = slau.jg[j] = i - 2;

j++;

slau.jg[j + 6] = slau.jg[j] = i - 1;

slau.jg[j + 7] = i;

i += 2;

j += 8;

}

}

void MKE::clearing()

{

int i;

for(i = 0; i < slau.n; i++)

slau.di[i] = slau.pr[i] = 0;

for(i = 0; i < slau.ig[slau.n]; i++)

slau.ggl[i] = slau.ggu[i] = 0;

}

void MKE::construction_global_SLAU()

{

clearing();

sigma *= omega;

Hi *= omega * omega;

for(int j = 0; j < ny - 1; j++)

for(int i = 0; i < nx - 1; i++)

{

definition_of_local_system(i, j);

}

conditions();

}

void MKE::definition_of_local_system(int tx, int ty)

{

double  hx = x[tx + 1] - x[tx],

hy = y[ty + 1] - y[ty],

hxy= hx / hy,

hyx= hy / hx,

hxhy = hx * hy;

pixel pix[4];

pix[0].x = x[tx  ]; pix[0].y = y[ty  ];

pix[1].x = x[tx+1]; pix[1].y = y[ty  ];

pix[2].x = x[tx  ]; pix[2].y = y[ty+1];

pix[3].x = x[tx+1]; pix[3].y = y[ty+1];

int i, j, k, l;

for(i = 0; i <8; i++) loc_v[i] = 0;

for(k = i = 0; i < 4; i++, k += 2)

for(l = j = 0; j < 4; j++, l += 2)

{

loc_m[k+1][l+1] = loc_m[k][l] = lambda * (B1[i][j] * hyx + B2[i][j] * hxy) - Hi * C[i][j] * hxhy;

loc_m[k+1][l] = sigma * C[i][j] * hxhy;

loc_m[k][l+1] = -loc_m[k+1][l];

loc_v[k]   += Fs(pix[j]) * C[i][j] * hxhy;

loc_v[k+1] += Fc(pix[j]) * C[i][j] * hxhy;

}