Современные вычислительные технологии

Второй этап - построение глобальной матрицы жёсткости и отыскание решения в перемещениях. Добавляем в проект новый класс - model с объявлением

class model
{
	// элементы
	std::vector<element*> _elements;
	// соответствие между узловым перемещениями элементов и индексов в глобальном векторе перемещений
	std::map<nodal_displacement, size_t> _nodes;
	// матрица жёсткости
	std::vector<std::vector<double>> _stiffness_matrix;
	// вспомогательные методы, для использования в конструкторе
	void create_dictionary();
	void create_the_stiffness_matrix();
public:
	/// <summary>
	/// конструктор
	/// </summary>
	/// <param name="elements">набор элементов</param>
	model(const std::vector<element*>& elements);

	/// <summary>
	/// отыскание решения в перемещениях
	/// </summary>
	/// <param name="constraints"> закрепления </param>
	/// <param name="loads"> нагрузки </param>
	/// <returns> решение задачи в перемещениях </returns>
	std::vector<double> solve(
		const std::map<nodal_displacement, double>& constraints,
		const std::map<nodal_displacement, double>& loads) const;
};

1. Прежде всего в классе необходимо создать открытый конструктор. Обратите внимание на поле _nodes. Это поле является ассоциативным контейнер std::map<nodal_displacement, size_t>. Для того, чтобы компилятор позволил вам создать такой ассоциативный контейнер, в классе nodal_displacement следует переопределить две операции сравнения следующим образом: 

           bool operator==(const nodal_displacement& that) const;
   	bool operator<(const nodal_displacement& that) const;

   (это перегрузка в виде метода класса, также при перегрузке возможно использование дружественных функций). При перегрузке помните о поле _number класса node, а также о том, что элементы перечислений можно сравнивать друг с другом непосредственно. Добавляйте открытые геттеры там, где их не хватает.
2. Поле _nodes задает соответствие между узловым перемещение и номером узлового перемещения в глобальном векторе узловых перемещений. Класс element содержит вектор, который возвращает вектор узловых перемещений. Он имеет столько же компонент, сколько степеней свободы имеет вектор (например, linear_element возвращает вектор из четырёх компонент). Для заполнения поля _nodes удобно создать вспомогательный метод create_dictionary(). Алгоритм его работы следующий:
   • заводим переменную типа size_t для счётчика
   • в цикле по элементам для каждого элемента запускаем цикл по узловым перемещениям.
   • если _nodes не содержит ключевое значение, соответствующее очередному узловому перемещению, добавляем в _nodes новую пару ключ-значение, при этом ключ равен очередному узловому перемещению, а значение - счётчику, сам счётчик при этом увеличивается на единицу.
3. Далее в конструкторе следует создать глобальную матрицу жёсткости. Размер матрицы жёсткости равен размеру _nodes. У каждого элемента есть своя матрица жёсткости в глобальной системе координат. Каждой строке (столбцу) этой матрицы соответствует своё узловое перемещение. Каждому из этих узловых перемещений соответствует номер, который может быть получен из таблицы _nodes. Таким образом, каждому индексу в локальной матрице соответствует индекс в глобальной матрице. Вспомогательный закрытый метод create_the_stiffness_matrix() должен делать следующее: 
• Создать матрицу жёсткости соответствующего размера и заполнить её нулями (при этом удобно воспользоваться конструктором вектора, принимающим размер)
• В цикле по элементам для каждого создать свою матрицу жёсткости и прибавить значения её коэффициентов к соответствующим коэффициентам глобальной матрицы жёсткости
4. Класс содержит указатели, следовательно, чтобы избежать утечек памяти, нужен нетривиальный деструктор. Деструктор в цикле проходит по вектору элементов и для каждого вызывает оператор delete. 
   Замечание. Если класс имеет нетривиальный деструктор, то он должен содержать также содержать явные копирующий конструктор и перегрузку оператора присвоения. Можно сделать их, а можно запретить копирование следующим образом:

   model(const model&) = delete;
   model& operator=(const model&) = delete;  

5. После того, как создана глобальная матрица жёсткости, создаём открытый метод

   std::vector<double> solve(
   	const std::map<nodal_displacement, double>& constraints,
   	const std::map<nodal_displacement, double>& loads) const;   

   Переменная constraints - ассоциативный контейнер, содержащий в качестве ключа узловые перемещения, а в качестве значения - значения узловых перемещений. Переменная loads содержит в качестве ключей узловые перемещения, а в качестве значений - нагрузки в направлении перемещений. Примерный алгоритм метода solve следующий:
• используя loads, создаём правую часть
• используя constraints, преобразуем правую часть. Предположим, что j - номер узлового перемещения, значение которого задано, i - индекс элемента правой части. Тогда элемент правой части изменяется следующим образом: \[ f_i=f_i-a_{ij}u_j, \] где \(f_i\) - элемент правой части, \(a_{ij}\) - элемент матрицы жёсткости, \(u_i\) - заданное значение узлового перемещения.
• решить СЛАУ \[ \sum\limits_{i=1}^na_{ij}u_j=f_i \] таким образом, чтобы в процессе решения не участвовали узловые перемещения, значения которых заданы (например, можно исключить из системы строки и столбцы матрицы, соответствующие узловым перемещениям, указанным в переменной constraints)
6. Осталось опробовать в деле то, что получилось. Напишем в main следующий код:

   int main() {
   	// создаём узлы
   	node n1 = { 1, 0.0, 3.0 };
   	node n2 = { 2, 0.0, 0.0 };
   	node n3 = { 3, 0.0, -3.0 };
   	node n4 = { 4, 4.0, 0.0 };
   	
           // создаём элементы
           element* el1 = new linear_element(n1, n4, 2e11, 0.1e-3);
   	element* el2 = new linear_element(n2, n4, 2e11, 0.1e-3);
   	element* el3 = new linear_element(n3, n4, 2e11, 0.1e-3);
   	model m = { {el1, el2, el3} };
   
   	// создаём закрепления
   	std::map<nodal_displacement, double> constraints =
   	{
   		{{n1, ux}, 0},
   		{{n1, uy}, 0},
   		{{n2, ux}, 0},
   		{{n2, uy}, 0},
   		{{n3, ux}, 0},
   		{{n3, uy}, 0},
   	};
   
   	// создаём нагрузки
   	std::map<nodal_displacement, double> loads =
   	{
   		{{n4, uy}, -1000}
   	};
   
   	// решаем СЛАУ
   	const auto solution = m.solve(constraints, loads);
   	for (auto x:solution)
   	{
   		std::cout << x << std::endl;
   	}
   
   	system("pause");
   }

   Если всё сделано правильно, вывод примерно такой (возможен другой порядок элементов вектора)

   0
   0
   0
   -0.000347222
   0
   0
   0
   0
   Для продолжения нажмите любую клавишу . . .