Расчет плоских ферм и арок. Расчетная схема фермы с принятой нумерацией узлов и стержней, страница 7

  }

//                                 

// Вывод результатов расчета в файл

//                                 

// Пункт 22

// Открыть файл "result.txt" на запись

  FILE *fileResult = fopen("result.txt", "w");

// Если файл не был открыт, то вывести сообщение об ошибке и выйти из

// программы

  if( fileResult == NULL )

  {

     cout << "Cannot open file" << endl;

     return 1;

  }

// Вывести в файл перемещения узлов

  fprintf(fileResult, "\nПеремещения узлов\n");

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

  {

     long place = 2*i;

     fprintf(fileResult, "%3d: x = %13.6lf ; y = %13.6lf\n", i,

Displace[ place ], Displace[ place+1 ]);

  }

// Вывести в файл усилия в элементах

  fprintf(fileResult, "\nУсилия в элементах\n");

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

  {

     long place = 2*i;

     fprintf(fileResult, "%3d: %15.3lf\n", i, rod[i].N);

  }

// Закрыть файл

  fclose(fileResult);

//                            

// Завершение работы программы

//                            

// Пункт 23

// Удаление выделенной памяти ...

// ... под объекты узлов, стержней, сил

  delete []node;

  delete []rod;

  delete []force;

// ... под глобальную матрицу жесткости

  for( i = 0; i<2*numNodes; i++ )

     delete []MatrixA[i];

  delete []MatrixA;

// ... под матрицу-столбец внешних сил

  delete []MatrixB;

// ... под матрицу-столбец перемещений узлов (с учетом связей)

  delete []MatrixX;

// ... под матрицу-столбец перемещений узлов

  delete []Displace;

// Выход из программы

  return 0;

}

Файл Solve.cpp - решение системы уравнений

#include "stdafx.h"

// Подключение библиотек математических функций, стандартного ввода-вывода и

// процессов

#include <math.h>

#include <stdio.h>

#include <process.h>

// Объявление функции вывода сообщений об ошибках при разложении системы

// уравнений (о нулевой или вырожденной матрице)

void Singular(int why, long iStr);

// Функция LU-разложения и решения разложенной системы уравнений. Суть

// переменных и кода этой функции подробно не комментируется, так как это

// выходит за рамки методического указания.

void LUDecomposition(long n, double **MatrixA, double *MatrixB, double *MatrixX)

{

// Разложение

  long *ps = new long[n];

  double *scales = new double[n];

  int i, j, k, pivotindex;

  double normrow, pivot, size, biggest, mult, ep=1e-12;

// Присвоение начальных значений массивам ps, LU, scales

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

     ps[i] = i;

     normrow = 0;

     for( j = 0; j < n; j ++ ) {

       if( normrow < fabs ( MatrixA[i][j] ) )

          normrow = fabs ( MatrixA[i][j] );

     }

     if( normrow != 0 )

       scales[i] = 1/normrow;

     else {

       scales[i] = 0;

       Singular(0, i);

     }

  }

// Метод исключения Гаусса с частичным упорядочиванием

  for( k = 0; k < n-1; k ++ ) {

     biggest = 0;

     for( i = k; i < n; i ++ ) {

       size = fabs( MatrixA[ps[i]][k] ) * scales[ps[i]];

       if( biggest < size ) {

          biggest = size;

          pivotindex = i;

       }

     }

     if( biggest == 0 ) {

       Singular(1, 0);

       goto endkloop;

     }

     if( pivotindex != k ) {

       j = ps[k];

       ps[k] = ps[pivotindex];

       ps[pivotindex] = j;

     }

     pivot = MatrixA[ps[k]][k];

    for( i = k+1; i < n; i ++ ) {

       MatrixA[ps[i]][k] = mult = MatrixA[ps[i]][k]/pivot;

       if( mult != 0 )

          for( j = k+1; j < n; j ++ )

            MatrixA[ps[i]][j] = MatrixA[ps[i]][j] - mult * MatrixA[ps[k]][j];

     }

     endkloop:;

  }

  if( MatrixA[ps[n-1]][n-1] < ep )

     Singular(1, 0);

  delete []scales;

// Решение разложенной системы уравнений

  double dot;

// Прямой ход

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

     dot = 0;

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

       dot += MatrixA[ps[i]][j] * MatrixX[j];

     MatrixX[i] = MatrixB[ps[i]] - dot;

  }

// Обратный ход

  for( i = n-1; i >= 0; i -- ) {

     dot = 0;

     for( j = i+1; j < n; j ++ )

       dot += MatrixA[ps[i]][j] * MatrixX[j];

     MatrixX[i] = ( MatrixX[i] - dot ) / MatrixA[ps[i]][i];

  }

  delete []ps;

}