Разработка программы решения гармонической задачи методом конечных элементов (Лабораторная работа № 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);

                     addition_in_SLAU(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;

              }