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

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

                     slau.di[k-1 + d] = slau.di[k-2 + d] = slau.di[k+1] = slau.di[k] += 2.0*Betta*hy / 6.;

                     slau.ggl[k] = slau.ggu[k] += Betta*hy / 6.;

                     }

              }

}

int MKE::save_information()

{

       FILE *file;

       int maxiter, i, j, k;

       double eps;

       pixel p;

       //     для LOS определяем max число итераций и eps относительной невязки

       if(!(file = fopen(Kuslau, "r")))return 1;

       fscanf(file, "%d %le", &maxiter, &eps);

       fclose(file);

       //     сохраняем данные для LOS в разреженно-строчном формате

       if(preservation_of_file(IG, 0, slau.ig, slau.n + 1))return 2;

       if(preservation_of_file(JG, 0, slau.jg, slau.ig[slau.n]))return 3;

       if(preservation_of_file(GGL, slau.ggl, 0, slau.ig[slau.n]))return 4;

       if(preservation_of_file(GGU, slau.ggu, 0, slau.ig[slau.n]))return 5;

       if(preservation_of_file(DI, slau.di, 0, slau.n))return 6;

       if(preservation_of_file(PR, slau.pr, 0, slau.n))return 7;

       if(!(file = fopen(KUSLAU, "w")))return 8;

       fprintf(file, "%d\n%d\n%le", slau.n, maxiter, eps);

       fclose(file);

       if(!(file = fopen(U, "w")))return 9;       //     сохраняем точное решение для сравнения с численным

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

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

              {

                     p.x = x[i];

                     p.y = y[j];

                     fprintf(file, "%.13le\n%.13le\n", Us(p), Uc(p));

              }

       fclose(file);

       //сохраняем данные для LU в профильном формате

       if(!(file = fopen(MATRIX, "w")))return 10;//      сохранение матрицы в профильном формате

       fprintf(file, "%d\n1 1 ", slau.n);

       for(k = 1, i = 1; i < slau.n; i++)//преобразование ig (индексация элементов начинается с 1) к пр. формату

       {

              j = i - slau.jg[slau.ig[i]]; //число элементов в строке

              k += j;

              fprintf(file, "%d ", k);

       }

       gg_gen_in_profile(file, slau.ggl);  //переводим матрицу (ggl, ggu) в проф. формат

       gg_gen_in_profile(file, slau.ggu);

       fprintf(file, "\n\n");

       for(i = 0; i < slau.n; i++)                //сохраняем диагональ

              fprintf(file, "%.13le ", slau.di[i]);

       fclose(file);

       if(preservation_of_file(VECTOR, slau.pr, 0, slau.n))return 11;

       if(!(file = fopen(XY, "w")))return 10;     //сохранение координат всех узлов сетки

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

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

                     fprintf(file, "%.10le\t%.10le\n%.10le\t%.10le\n", x[j], y[i], x[j], y[i]);

       fclose(file);

       return 0;

}

void MKE::gg_gen_in_profile(FILE *file, double *gg)

{

       int i, j, k;

       fprintf(file, "\n");

       for(i = 1; i < slau.n; i++) //цикл по строкам матрицы

       {

              fprintf(file, "\n");

              for(j = slau.ig[i]; j < slau.ig[i+1] - 1; j++)    //цикл по элементам строки

              {

                     fprintf(file, "%.13le ", gg[j]);

                     for(k = 1; slau.jg[j] + k < slau.jg[j+1]; k++)       //определяем сколько не хватает нулей для перевода

                            fprintf(file, "0. ");

              }

              fprintf(file, "%.13le ", gg[j]);    //последний элемент в строке

              for(k = 1; slau.jg[j] + k < i; k++)

                     fprintf(file, "0. ");

       }

}

int MKE::preservation_of_file(char *s, double *le, int *d, int n)

{

       FILE *file;

       if(!(file = fopen(s, "w")))return 1;

       if(le == NULL) //     определяем целые или вещественные данные

              for(int i = 0; i < n; i++)

                     fprintf(file, "%d ", d[i]);