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

}

void MKE::addition_in_SLAU(int ix, int jy)

{

       int i, j, k,

              n[8];  //     номера вершин локального элемента

       //     определяем эти номера

       n[0] = (jy * nx + ix) * 2;

       n[4] = n[0] + nx * 2;

       for(i = 1; i < 4; i++)

       {

              n[i] = n[0] + i;

              n[i + 4] = n[4] + i;

       }

       //     добавление элементов локальной системы в глобальную

       for(i = 0; i < 8; i++)

       {

              slau.di[n[i]] += loc_m[i][i];

              slau.pr[n[i]] += loc_v[i];

              for(j = 0; j < i; j++)

              {

                     k = slau.ig[n[i]];   //     указывает на первый элемент в строке

                     while(slau.jg[k] != n[j]) k++;      //     поиск места элемента в массиве gg

                     slau.ggl[k] += loc_m[i][j];

                     slau.ggu[k] += loc_m[j][i];

              }

       }

}

//     учитываем краевые условия

void MKE::conditions()

{

       int i, k, d = (ny - 1) * nx * 2;

       //     первые краевые сверху и снизу

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

       {

              slau.di[k+1 + d] = slau.di[k + d] = slau.di[k+1] = slau.di[k] = BIG;

              slau.pr[k]   = BIG * C1sYBottom(x[i]);

              slau.pr[k+1] = BIG * C1cYBottom(x[i]);

              slau.pr[k + d]   = BIG * C1sYTop(x[i]);

              slau.pr[k+1 + d] = BIG * C1cYTop(x[i]);

       }

       d = 2 * nx;

       if(CN() == 1)

              //     первые краевые справа и слева

              for(k = i = 0; i < ny; i++, k += d)

              {

                     slau.di[k-1 + d] = slau.di[k-2 + d] = slau.di[k+1] = slau.di[k] = BIG;

                     slau.pr[k]   = BIG * C1sXLeft(y[i]);

                     slau.pr[k+1] = BIG * C1cXLeft(y[i]);

                     slau.pr[k-2 + d] = BIG * C1sXRight(y[i]);

                     slau.pr[k-1 + d] = BIG * C1cXRight(y[i]);

              }

       else

              //     вторые краевые справа и слева

              if(CN() == 2)

              {

                     for(k = i = 0; i < ny - 1; i++, k += d)

                     {

                            double hy = y[i+1] - y[i];

                            slau.pr[k]     += (2.* C2sXLeft(y[i+1]) + C2sXLeft(y[i])) * hy / 6.;

                            slau.pr[k+1]   += (2.* C2cXLeft(y[i+1]) + C2cXLeft(y[i])) * hy / 6.;

                            slau.pr[k+d]   += (C2sXLeft(y[i]) + 2.* C2sXLeft(y[i+1])) * hy / 6.;

                            slau.pr[k+d+1] += (C2cXLeft(y[i]) + 2.* C2cXLeft(y[i+1])) * hy / 6.;

                            slau.pr[k-2+d] += (C2sXRight(y[i]) + 2.* C2sXRight(y[i+1])) * hy / 6.;

                            slau.pr[k-1+d] += (C2cXRight(y[i]) + 2.* C2cXRight(y[i+1])) * hy / 6.;

                            slau.pr[k-2+2*d] += (2.* C2sXRight(y[i]) + C2sXRight(y[i+1])) * hy / 6.;

                            slau.pr[k-1+2*d] += (2.* C2cXRight(y[i]) + C2cXRight(y[i+1])) * hy / 6.;

                     }

              }

              //     третьи краевые справа и слева

              else

              {

                     for(k = i = 0; i < ny - 1; i++, k += d)

                     {

                            double hy = y[i+1] - y[i];

                            slau.pr[k]     += (2.* C3sXLeft(y[i+1]) + C3sXLeft(y[i])) *Betta*hy / 6.;

                            slau.pr[k+1]   += (2.* C3cXLeft(y[i+1]) + C3cXLeft(y[i])) *Betta*hy / 6.;

                            slau.pr[k+d]   += (C3sXLeft(y[i]) + 2.* C3sXLeft(y[i+1])) *Betta*hy / 6.;

                            slau.pr[k+d+1] += (C2cXLeft(y[i]) + 2.* C3cXLeft(y[i+1])) *Betta*hy / 6.;

                            slau.pr[k-2+d] += (C3sXRight(y[i]) + 2.* C3sXRight(y[i+1])) *Betta*hy / 6.;

                            slau.pr[k-1+d] += (C3cXRight(y[i]) + 2.* C3cXRight(y[i+1])) *Betta*hy / 6.;

                            slau.pr[k-2+2*d] += (2.* C3sXRight(y[i]) + C3sXRight(y[i+1])) *Betta*hy / 6.;