[румата]

Эта тема создана с целью поиска и сбора информации по любым из существующих способов и методам определения расчетных длин, а также с целью их обсуждения, критики или формулировок своих собственный методик применительно к расчетам на устойчивость и проверке по предельной гибкости согласно норм РФ и, возможно, других стран.

[Бам]

Я так понимаю, по всей стране всё расчеты сжатых рам нужно немедленно прервать, пока ваши исследования не будут закончены? Нас ждут сенсационные результаты?

[Ильнур]

Сообщение от Бахил ... А длины то будут?

Это несходящийся процесс.

Нас ждут сенсационные результаты?

Да, будет показано, что правильная расчетная длина - величина неопределяемая в принципе, если стержней более 1 шт.

[Бахил]

Сообщение от Ильнур правильная расчетная длина - величина неопределяемая в принципе, если стержней более 1 шт.

Да! Это надо отлить в бронзе и поместить на фронтоны всех строительных вузов.

[crossing]

Надо вот этого вот румату забанить на всегда. Как забанили этого семсонофф и онже илюха. И румата из тогоже набора.

[румата]

Сообщение от crossing Надо вот этого вот румату забанить на всегда.

Offtop: Интересно, а за что надо банить румату?

[Бам]

Offtop:

Сообщение от румата Offtop: Интересно, а за что надо банить румату?

"А вы не любите пролетариат..."(с)

[crossing]

Сообщение от румата Offtop: Интересно, а за что надо банить румату?

За то, что не понимаешь о чём спрашиваешь. Вот ты тему создал, а у тебя в мозгах контекста нет.
Ты через полгода ещё такуюже сосдаш.

[румата]

Сообщение от crossing За то, что не понимаешь о чём спрашиваешь.

В теме вообще нет вопроса. Эта тема для сбора информации и поиска способов и методов определния расчетных длин пригодных для расчетов по СП16.13330.

Сообщение от crossing Вот ты тему создал, а у тебя в мозгах контекста нет.

Снова интересно. Что ты подразумеваешь под контекстом, которого у меня нет в мозгах?

Сообщение от crossing Ты через полгода ещё такуюже сосдаш.

Зачем? Мне этой вполне достатотчно на всю оставшуюся жизнь.

[Нубий-IV]

То, что страшный геомнелин - на самом деле обычная сумма перемещений с давно известным коэффициентом 1/(1-N/Ncr), немного упрощает задачу. Были бы формы независимы - можно было бы просто суммировать всех, кто дает положительный вклад, как при линиях влияния. Проблема в том, что формы конфликтуют между собой за право вписаться в допуск. И возможна ситуация, когда главную по влиянию форму можно еще увеличить, добавив немного отрицательной второй формы.

Например, для стержня Эйлера при учете пары первых форм:

• В первом узле коэффициенты формы связаны условием вписаться в допуск:
  
  
• Аналогично, во втором узле:
  
  
• Если перебрать все узлы - получится область допустимых значений.
  
  
• Если на эту область нанести изолинии постоянных N, можно найти наихудшую комбинацию:
  
  
  Сейчас изолинии почти вертикальные, но это особенность одиночного стержня, где вторая критическая в разы выше первой. В рамах, если соседние формы почти равны, линии наклонятся ближе к 45°; а если и соотношение форм удачно ляжет - то наклонятся еще сильнее. Линии не параллельны, они постепенно выходят на вертикаль при приближении к Ncr1, или к фиксированному углу, если RA<Ncr1.
• 
  Сразу можно отбраковать только случай двух отрицательных коэффициентов (он означает максимальное напряжение в противоположной точке сечения).
  В первой четверти решение есть всегда. В четвертях с разными знаками возможны области, где решения нет. Но в целом при других соотношениях форм и критических сил возможны и варианты, когда надо выбирать единственную форму, и варианты, когда хуже оказываются комбинации. Например, чудо-рама - случай выбора второй формы, а Эйлер с постоянным допуском - случай выбора комбинации.
• С учетом того, что зависимость N от коэффициентов формы нелинейная, получается задача нелинейного программирования. Сходу можно только сказать, что минимум попадает в узел области, потому что изолинии N прямые, то есть достаточно перебора конечного числа точек.
• Когда коэффициенты форм выбраны, поиск N - это поиск наименьшего корня полинома. Каждая форма добавляет единицу к степени уравнения. При одной форме уравнение квадратное, при двух - кубическое и т.д.
• Кто может сказать - проще это, чем поиск минимума среди всех возможных искривлений через последовательное решение системы уравнений MKЭ?
• В случае с Эйлером подсказку дает перемножение эпюр, это намного проще, чем искать минимум бесконечного числа форм. Для рамы возможны какие-нибудь аналогичные соображения?

Сообщение от Бам по всей стране всё расчеты сжатых рам нужно немедленно прервать, пока ваши исследования не будут закончены?

Риторические вопросы слишком сложны, их надо постить в раздел форума "философия". А продемонстрировать острый ум можно намного проще. Например, предоставить готовый ответ на вопрос темы. Или хотя бы ткнуть в какую-нибудь ошибку (в посте 457 есть минимум одна; не буду удалять, пусть лежит).

Сообщение от Ильнур расчетная длина - величина неопределяемая в принципе, если стержней более 1 шт

Сообщение от Бахил Да! Это надо отлить в бронзе и поместить на фронтоны всех строительных вузов.

Сообщение от румата crossing Ты через полгода ещё такуюже сосдаш. румата Зачем?

Если у вас нету темы,
В нее не нагадит сосед.
И не обсуждать вам проблемы
С тем, для кого,
С тем, для кого,
С тем, для кого их нет.

[Бахил]

"Формы" линейно независимы и ортогональны. Фундаментальное свойство собственных чисел и векторов.

[Нубий-IV]

Сообщение от Бахил "Формы" линейно независимы и ортогональны.

Это в задаче расчета на устойчивость так: сумма сил в узле на дополнительных отклонениях равна нулю, откуда ненулевое решение при нулевом векторе нагрузок - это только собственный вектор.

А в геомнелине сюда добавляется начальное искривление и кривые нагрузки, которые добавляют усилий на дополнительных отклонениях и дают справа ненулевой вектор. В такой задаче вектор неизвестных надо разбирать на собственные, как в сейсмике или пульсациях. И находить коэффициенты формы для каждого вектора - из условия, чтобы в сумме заданное начальное искривление получилось. А начальное искривление, очевидно, должно в допуск вписаться.

Затея была сверить СП с честным нелинейным расчетом. Одна мелочь мешает - в геомнелин надо задать правильное начальное искривление. Правильное - это при котором несущая будет минимальной. Пример чудо-рамы показывает, что просто взять первую форму недостаточно - можно пролететь на 10-15%. Анализ изополей фи из СП показывает, что даже выбор между двумя формами может дать до 25% разницы, т.е. чудо-рама еще не худший вариант. А пример со стержнем Эйлера показывает, что и форм не всегда двух хватит. Как получится из геомнелина минимум доставать - можно будет сверять расчетные длины.

Пока что вижу только решение "в лоб" - свалять свой микро-мкэ для плоских рам, и в нем поиск минимума по начальным искривлениям организовать для схем из СП. Но запрограммить это небыстро будет.

[Нубий-IV]

Заготовка для обработки напильником. Исходник C++ (собирать как консольную программу в VS или CodeBlocks):

Код:

[Выделить все]

 
#include <string>
#include <sstream>
#include <iostream>
#include <fstream>
#include <list>
#include <math.h>
#include <assert.h>

using namespace std;

// Число степеней свободы в узле: X = 0, Y = 1, Rz = 2
const size_t NODE_NUM_DOFS = 3;
// Число узлов конечного элемента - двухузловой стержень
const size_t ELEMENT_NUM_NODES = 2;
// Число сечений по длине элемента для вычисления усилий и напряжений
const size_t ELEMENT_NUM_SECTIONS = 5;
// Размер матрицы жесткости конечного элемента
const size_t ELEMENT_MATRIX_SIZE = ELEMENT_NUM_NODES * NODE_NUM_DOFS;

// Простейший массив для хранения элементов модели, матриц жесткости и результатов
template <typename T> class Array {
public:
	Array(): size(0), value(nullptr) {}

	// Выделение памяти
	void Allocate(size_t size)
	{
		assert(size > 0);
		assert(this->size == 0);
		assert(this->value == nullptr);

		this-> size = size;
		this->value = new T[size];
	}

	size_t Size()
	{
		return size;
	}

	// Оператор доступа по индексу
	T& operator()(size_t i)
	{
		assert(size != 0);
		assert(value != nullptr);
		assert(i < size);

		return value[i];
	}
private:
	size_t size;
	T* value;
};

// Узел в модели
struct Node {
	double x;
	double y;
};

// Материал в модели
struct Material {
	double A; // Площадь
	double W; // Момент сопротивления
	double I; // Момент инерции
	double E; // Модуль упругости
};

// Элемент в модели
struct Element {
	// Номер материала
	size_t iMaterial;
	// Номера узлов
	size_t iNode[ELEMENT_NUM_NODES];
};

// Связь в модели
struct Support {
	// Номер узла
	size_t iNode;
	// Номер направления связи
	size_t iDirection;
};

// Нагрузка в модели
struct Force {
	size_t iNode;
	size_t iDirection;
	double value;
};

// Матрица жесткости элемента
class ElementMatrix {
public:
	// При создании - обнуление
	ElementMatrix()
	{
		SetZero();
	}

	// Оператор доступа по индексам
	double& operator()(size_t i, size_t j)
	{
		assert(i < ELEMENT_MATRIX_SIZE);
		assert(j < ELEMENT_MATRIX_SIZE);

		return value[i][j];
	}

	// Возвращает локальный номер узла по индексу элемента матрицы
	size_t LocalNodeId(size_t i)
	{
		assert(i < ELEMENT_MATRIX_SIZE);

		return i / NODE_NUM_DOFS;
	}

	// Возвращает номер направления по индексу элемента матрицы
	size_t DirectionId(size_t i)
	{
		assert(i < ELEMENT_MATRIX_SIZE);

		return i % NODE_NUM_DOFS;
	}

	// Обнуление всех элементов
	void SetZero()
	{
		for(size_t i = 0; i < ELEMENT_MATRIX_SIZE; i++)
			for(size_t j = 0; j < ELEMENT_MATRIX_SIZE; j++)
				value[i][j] = 0.0;
	}
private:
	double value[ELEMENT_MATRIX_SIZE][ELEMENT_MATRIX_SIZE];
};

// Матрица поворотов элемента
class TransformationMatrix {
public:
	// При создании - обнуление
	TransformationMatrix()
	{
		SetZero();
	}

	// Оператор доступа по индексам
	double& operator()(size_t i, size_t j)
	{
		assert(i < ELEMENT_MATRIX_SIZE);
		assert(j < ELEMENT_MATRIX_SIZE);

		return value[i][j];
	}

	// Обнуление всех элементов
	void SetZero()
	{
		for(size_t i = 0; i < ELEMENT_MATRIX_SIZE; i++)
			for(size_t j = 0; j < ELEMENT_MATRIX_SIZE; j++)
				value[i][j] = 0;
	}
private:
	double value[ELEMENT_MATRIX_SIZE][ELEMENT_MATRIX_SIZE];
};

// Вектор элемента - перемещения, реакции в узлах
class ElementVector {
public:
	// При создании - обнуление
	ElementVector()
	{
		SetZero();
	}

	// Оператор доступа по индексу
	double& operator()(size_t i)
	{
		assert(i < ELEMENT_MATRIX_SIZE);

		return value[i];
	}

	// Возвращает локальный номер узла по индексу элемента вектора
	size_t LocalNodeId(size_t i)
	{
		assert(i < ELEMENT_MATRIX_SIZE);

		return i / NODE_NUM_DOFS;
	}

	// Возвращает номер направления по индексу элемента вектора
	size_t DirectionId(size_t i)
	{
		assert(i < ELEMENT_MATRIX_SIZE);

		return i % NODE_NUM_DOFS;
	}

	// Обнуление всех элементов
	void SetZero()
	{
		for(size_t i = 0; i < ELEMENT_MATRIX_SIZE; i++)
			value[i] = 0.0;
	}
private:
	double value[ELEMENT_MATRIX_SIZE];
};

// Вектор элемента - усилия и напряжения в расчетных сечениях
class ElementResult {
public:
	// При создании - обнуление
	ElementResult()
	{
		SetZero();
	}

	// Оператор доступа по индексу
	double& operator()(size_t i)
	{
		assert(i < ELEMENT_NUM_SECTIONS);

		return value[i];
	}

	// Обнуление всех элементов
	void SetZero()
	{
		for(size_t i = 0; i < ELEMENT_NUM_SECTIONS; i++)
			value[i] = 0.0;
	}
private:
	double value[ELEMENT_NUM_SECTIONS];
};

// Мартица коэффициентов системы - матрица жесткости
class SystemMatrix {
public:
	SystemMatrix(): size(0), value(nullptr) {}

	// Выделение памяти и обнуление
	void Allocate(size_t size)
	{
		assert(size > 0);
		assert(this->size == 0);
		assert(this->value == nullptr);

		this->size = size;
		this->value = new double*[size];
		for(size_t i = 0; i < size; i++)
			this->value[i] = new double[size];

		SetZero();
	}

	size_t Size()
	{
		return size;
	}

	// Оператор доступа по номеру узла и номеру перемещения
	double& operator()(size_t iNode, size_t iDirection, size_t jNode, size_t jDirection)
	{
		assert(size != 0);
		assert(value != nullptr);
		assert(iNode * NODE_NUM_DOFS + iDirection < size);
		assert(jNode * NODE_NUM_DOFS + jDirection < size);

		return value[iNode * NODE_NUM_DOFS + iDirection][jNode * NODE_NUM_DOFS + jDirection];
	}

	// Оператор доступа по индексу
	double& operator()(size_t i, size_t j)
	{
		assert(size != 0);
		assert(value != nullptr);
		assert(i < size);
		assert(j < size);

		return value[i][j];
	}

	// Обнуление всех элементов
	void SetZero()
	{
		for(size_t i = 0; i < size; i++)
			for(size_t j = 0; j < size; j++)
				this->value[i][j] = 0.0;
	}
private:
	size_t size;
	double** value;
};

// Вектор системы - нагрузки и неизвестные
class SystemVector {
public:
	SystemVector(): size(0), value(nullptr) {}

	// Выделение памяти и обнуление
	void Allocate(size_t size)
	{
		assert(size > 0);
		assert(this->size == 0);
		assert(this->value == nullptr);

		this-> size = size;
		this->value = new double[size];

		SetZero();
	}

	size_t Size()
	{
		return size;
	}

	// Оператор доступа по номеру узла и номеру перемещения
	double& operator()(size_t iNode, size_t iDirection)
	{
		assert(size != 0);
		assert(value != nullptr);
		assert(iNode * NODE_NUM_DOFS + iDirection < size);

		return value[iNode * NODE_NUM_DOFS + iDirection];
	}

	// Оператор доступа по индексу
	double& operator()(size_t i)
	{
		assert(size != 0);
		assert(value != nullptr);
		assert(i < size);

		return value[i];
	}

	// Обнуление всех элементов
	void SetZero()
	{
		for(size_t i = 0; i < size; i++)
			this->value[i] = 0.0;
	}
private:
	size_t size;
	double* value;
};

// Модель, загружаемая из файла
class Model {
public:
	Array<Node>     Nodes;
	Array<Material> Materials;
	Array<Element>  Elements;
	Array<Support>  Supports;
	Array<Force>    Forces;
};

// Данные конечных элементов для модели
class ElementsData {
public:
	// Матрицы жесткости в локальных осях
	Array<ElementMatrix> LocalMatrices;
	// Матрицы поворота
	Array<TransformationMatrix> TransformationMatrices;
	// Матрицы жесткости в глобальных осях
	Array<ElementMatrix> GlobalMatrices;

	// Перемещения элементов в локальных осях
	Array<ElementVector> LocalDofs;
	// Перемещения элементов в глобальных осях
	Array<ElementVector> GlobalDofs;

	// Реакции в узлах элементов в локальных осях
	Array<ElementVector> LocalReactions;

	// Продольные силы в расчетных сечениях элементов
	Array<ElementResult> LongitudinalForces;
	// Поперечные силы в расчетных сечениях элементов
	Array<ElementResult> ShearForces;
	// Моменты в расчетных сечениях элементов
	Array<ElementResult> BendingMoments;
	// Напряжения в расчетных сечениях элементов
	Array<ElementResult> Stresses;
public:
	// При создании - выделени памяти для данных
	ElementsData(Model& model): model(model)
	{
		LocalMatrices.Allocate (model.Elements.Size());
		TransformationMatrices.Allocate(model.Elements.Size());
		GlobalMatrices.Allocate(model.Elements.Size());

		LocalDofs.Allocate(model.Elements.Size());
		GlobalDofs.Allocate(model.Elements.Size());

		LocalReactions.Allocate(model.Elements.Size());

		LongitudinalForces.Allocate(model.Elements.Size());
		ShearForces.Allocate(model.Elements.Size());
		BendingMoments.Allocate(model.Elements.Size());
		Stresses.Allocate(model.Elements.Size());
	}

	// Вычисление матриц жесткости всех КЭ в локальных осях
	void CalcLocalMatrices()
	{
		for(size_t i = 0; i < LocalMatrices.Size(); i++)
			CalcLocalMatrix(i, LocalMatrices(i));
	}

	// Вычисление матриц поворота всех КЭ
	void CalcTransformationMatrices()
	{
		for(size_t i = 0; i < TransformationMatrices.Size(); i++)
			CalcTransformationMatrix(i, TransformationMatrices(i));
	}

	// Перевод матриц жесткости всех КЭ в глобальные оси осях
	void CalcGlobalMatrices()
	{
		for(size_t i = 0; i < GlobalMatrices.Size(); i++)
			CalcGlobalMatrix(LocalMatrices(i), TransformationMatrices(i), GlobalMatrices(i));
	}

	// Перевод перемещений всех КЭ в локальные оси
	void CalcLocalDofs()
	{
		for(size_t i = 0; i < LocalDofs.Size(); i++)
			CalcLocalDofs(GlobalDofs(i), TransformationMatrices(i), LocalDofs(i));
	}

	// Вычисление узловых реакций всех КЭ
	void CalcLocalReactions()
	{
		for(size_t i = 0; i < LocalReactions.Size(); i++)
			CalcLocalReactions(LocalMatrices(i), LocalDofs(i), LocalReactions(i));
	}

	// Вычисление продольных сил в расчетных сечениях всех КЭ
	void CalcLongitudinalForces()
	{
		for(size_t i = 0; i < LocalReactions.Size(); i++)
			CalcLongitudinalForces(LocalReactions(i), LongitudinalForces(i));
	}

	// Вычисление поперечных сил в расчетных сечениях всех КЭ
	void CalcShearForces()
	{
		for(size_t i = 0; i < ShearForces.Size(); i++)
			CalcShearForces(LocalReactions(i), ShearForces(i));
	}

	// Вычисление моментов сил в расчетных сечениях всех КЭ
	void CalcBendingMoments()
	{
		for(size_t i = 0; i < BendingMoments.Size(); i++)
			CalcBendingMoments(i, LocalReactions(i), BendingMoments(i));
	}

	// Вычисление напряжений в расчетных сечениях всех КЭ
	void CalcStresses()
	{
		for(size_t i = 0; i < Stresses.Size(); i++)
			CalcStresses(model.Materials(i), LongitudinalForces(i), ShearForces(i), BendingMoments(i), Stresses(i));
	}
private:
	Model& model;
private:
	// Вычисление матрицы жесткости КЭ
	void CalcLocalMatrix(size_t iElement, ElementMatrix& matrix)
	{
		Element element = model.Elements(iElement);
		Material material = model.Materials(element.iMaterial);
		Node node0 = model.Nodes(element.iNode[0]);
		Node node1 = model.Nodes(element.iNode[1]);

		double EA = material.E * material.A;
		double EI = material.E * material.I;

		double dx = node1.x - node0.x;
		double dy = node1.y - node0.y;

		double L1  = sqrt(dx * dx + dy * dy);
		double L2 = L1 * L1;
		double L3 = L2 * L1;

		double v____EA_L1 =      EA / L1;
		double v_12_EI_L3 = 12 * EI / L3;
		double v__6_EI_L2 =  6 * EI / L2;
		double v__4_EI_L1 =  4 * EI / L1;
		double v__2_EI_L1 =  2 * EI / L1;

		matrix(0,0) =   v____EA_L1;
		matrix(0,3) = - v____EA_L1;

		matrix(1,1) =   v_12_EI_L3;
		matrix(1,2) =   v__6_EI_L2;
		matrix(1,4) = - v_12_EI_L3;
		matrix(1,5) =   v__6_EI_L2;

		matrix(2,1) =   v__6_EI_L2;
		matrix(2,2) =   v__4_EI_L1;
		matrix(2,4) = - v__6_EI_L2;
		matrix(2,5) =   v__2_EI_L1;

		matrix(3,0) = - v____EA_L1;
		matrix(3,3) =   v____EA_L1;

		matrix(4,1) = - v_12_EI_L3;
		matrix(4,2) = - v__6_EI_L2;
		matrix(4,4) =   v_12_EI_L3;
		matrix(4,5) = - v__6_EI_L2;

		matrix(5,1) =   v__6_EI_L2;
		matrix(5,2) =   v__2_EI_L1;
		matrix(5,4) = - v__6_EI_L2;
		matrix(5,5) =   v__4_EI_L1;
	}

	// Вычисление матрицы поворота КЭ
	void CalcTransformationMatrix(size_t iElement, TransformationMatrix& matrix)
	{
		Element element = model.Elements(iElement);
		Node node0 = model.Nodes(element.iNode[0]);
		Node node1 = model.Nodes(element.iNode[1]);

		double dx = node1.x - node0.x;
		double dy = node1.y - node0.y;
		double L  = sqrt(dx * dx + dy * dy);

		double vSin = dy / L;
		double vCos = dx / L;

		matrix(0,0) =   vCos;
		matrix(0,1) =   vSin;

		matrix(1,0) = - vSin;
		matrix(1,1) =   vCos;

		matrix(2,2) =   1.0;

		matrix(3,3) =   vCos;
		matrix(3,4) =   vSin;

		matrix(4,3) = - vSin;
		matrix(4,4) =   vCos;

		matrix(5,5) =   1.0;
	}

	// Перевод матрицы жесткости КЭ
	void CalcGlobalMatrix(ElementMatrix& localMatrix, TransformationMatrix& tMatrix, ElementMatrix& globalMatrix)
	{
		static double tmp[ELEMENT_MATRIX_SIZE][ELEMENT_MATRIX_SIZE];

		for(size_t i = 0; i < ELEMENT_MATRIX_SIZE; i++)
			for(size_t j = 0; j < ELEMENT_MATRIX_SIZE; j++) {
				tmp[i][j] = 0.0;
				for(size_t k = 0; k < ELEMENT_MATRIX_SIZE; k++)
					tmp[i][j] += localMatrix(i,k) * tMatrix(k,j);
			}

		for(size_t i = 0; i < ELEMENT_MATRIX_SIZE; i++)
			for(size_t j = 0; j < ELEMENT_MATRIX_SIZE; j++) {
				globalMatrix(i,j) = 0.0;
				for(size_t k = 0; k < ELEMENT_MATRIX_SIZE; k++)
					globalMatrix(i,j) += tMatrix(k,i) * tmp[k][j];
			}
	}

	// Перевод перемещений КЭ в локальные оси
	void CalcLocalDofs(ElementVector& globalDofs, TransformationMatrix& tMatrix, ElementVector& localDofs)
	{
		for(size_t i = 0; i < ELEMENT_MATRIX_SIZE; i++) {
			localDofs(i) = 0.0;
			for(size_t k = 0; k < ELEMENT_MATRIX_SIZE; k++)
				localDofs(i) += tMatrix(i,k) * globalDofs(k);
		}
	}

	// Вычисление узловых реакций КЭ
	void CalcLocalReactions(ElementMatrix& localMatrix, ElementVector& localDofs, ElementVector& localReactions)
	{
		for(size_t i = 0; i < ELEMENT_MATRIX_SIZE; i++) {
			localReactions(i) = 0.0;
			for(size_t k = 0; k < ELEMENT_MATRIX_SIZE; k++)
				localReactions(i) += localMatrix(i,k) * localDofs(k);
		}
	}

	// Вычисление продольных сил в расчетных сечениях КЭ
	void CalcLongitudinalForces(ElementVector& localReactions, ElementResult& longitudinalForces)
	{
		for(size_t i = 0; i < ELEMENT_NUM_SECTIONS; i++)
			longitudinalForces(i) = -localReactions(0);
	}

	// Вычисление поперечных сил в расчетных сечениях КЭ
	void CalcShearForces(ElementVector& localReactions, ElementResult& shearForces)
	{
		for(size_t i = 0; i < ELEMENT_NUM_SECTIONS; i++)
			shearForces(i) = localReactions(1);
	}

	// Вычисление моментов сил в расчетных сечениях КЭ
	void CalcBendingMoments(size_t iElement, ElementVector& localReactions, ElementResult& bendingMoments)
	{
		Element element = model.Elements(iElement);

		Node node0 = model.Nodes(element.iNode[0]);
		Node node1 = model.Nodes(element.iNode[1]);

		double dx = node1.x - node0.x;
		double dy = node1.y - node0.y;

		double L = sqrt(dx * dx + dy * dy);

		for(size_t i = 0; i < ELEMENT_NUM_SECTIONS; i++) {
			double k = double(i) / double(ELEMENT_NUM_SECTIONS - 1);

			bendingMoments(i) = k * L * localReactions(1) - localReactions(2);
		}
	}

	// Вычисление напряжений в расчетных сечениях КЭ
	void CalcStresses(Material& material, ElementResult& longitudinalForces, ElementResult& shearForces, ElementResult& bendingMoments, ElementResult& stresses)
	{
		for(size_t i = 0; i < ELEMENT_NUM_SECTIONS; i++)
			stresses(i) = longitudinalForces(i) / material.A + bendingMoments(i) / material.W;
	}
};

// Решатель. Матрица жесткости системы, вектора нагрузок и неизвестных
class Solver {
public:
	// Матрица жесткости системы
	SystemMatrix SystemMatrix;
	// Вектор нагрузок системы
	SystemVector SystemForces;
	// Вектор неизвестных
	SystemVector SystemDofs;
public:
	// При создании - выделение памяти
	Solver(Model& model, ElementsData& elementsData):
		model(model), elementsData(elementsData)
	{
		SystemMatrix.Allocate(model.Nodes.Size() * NODE_NUM_DOFS);
		SystemForces.Allocate(model.Nodes.Size() * NODE_NUM_DOFS);
		SystemDofs.Allocate(model.Nodes.Size() * NODE_NUM_DOFS);
	}

	// Сборка матрицы жесткости системы
	void PlaceElementGlobalMatrices()
	{
		for(size_t i = 0; i < model.Elements.Size(); i++)
			PlaceElementGlobalMatrix(i, elementsData.GlobalMatrices(i));
	}

	// Установка нагрузок
	void PlaceForces()
	{
		for(size_t i = 0; i < model.Forces.Size(); i++)
			PlaceForce(model.Forces(i));
	}

	// Установка связей
	void PlaceSupports()
	{
		for(size_t i = 0; i < model.Supports.Size(); i++)
			PlaceSupport(model.Supports(i));
	}

	// Решение системы уравнений
	void SolveGauss()
	{
		for(size_t i = 0; i < SystemForces.Size(); i++)
			SystemDofs(i) = SystemForces(i);

		// Прямой ход
		for(size_t i = 1; i < SystemMatrix.Size(); i++)
			for(size_t n = 0; n < i; n++) {
				double k = SystemMatrix(n,i) / SystemMatrix(n,n);
				//for(size_t j = i; j < SystemMatrix.Size(); j++)
				for(size_t j = 0; j < SystemMatrix.Size(); j++)
					SystemMatrix(i,j) -= k * SystemMatrix(n,j);
				SystemDofs(i) -= k * SystemDofs(n);
			}

		// Обратный ход
		for(size_t i = SystemMatrix.Size()-1; i < SystemMatrix.Size(); i--) {
			for(size_t j = i + 1; j < SystemMatrix.Size(); j++)
				SystemDofs(i) -= SystemMatrix(i,j) * SystemDofs(j);
			SystemDofs(i) /= SystemMatrix(i,i);
		}
	}

	// Выгрузка перемещений системы в вектора перемещений КЭ
	void PlaceElementGlobalDofs()
	{
		for(size_t i = 0; i < model.Elements.Size(); i++)
			PlaceElementGlobalDofs(i, elementsData.GlobalDofs(i));
	}
private:
	Model& model;
	ElementsData& elementsData;
private:
	// Установка матрицы жесткости КЭ в матрицу системы
	void PlaceElementGlobalMatrix(size_t iElement, ElementMatrix& elementGlobalMatrix)
	{
		Element element = model.Elements(iElement);

		for(size_t i = 0; i < ELEMENT_MATRIX_SIZE; i++)
			for(size_t j = 0; j < ELEMENT_MATRIX_SIZE; j++) {
				size_t iLocalNode = elementGlobalMatrix.LocalNodeId(i);
				size_t jLocalNode = elementGlobalMatrix.LocalNodeId(j);

				size_t iGlobalNode = element.iNode[iLocalNode];
				size_t jGlobalNode = element.iNode[jLocalNode];

				size_t iDirection = elementGlobalMatrix.DirectionId(i);
				size_t jDirection = elementGlobalMatrix.DirectionId(j);

				SystemMatrix(iGlobalNode, iDirection, jGlobalNode, jDirection) += elementGlobalMatrix(i,j);
			}

	}

	// Установка нагрузки в вектор системы
	void PlaceForce(Force& force)
	{
		size_t I = force.iNode * NODE_NUM_DOFS + force.iDirection;
		SystemForces(I) = force.value;
	}

	// Установка связи  в матрицу системы
	void PlaceSupport(Support& support)
	{
		size_t I = support.iNode * NODE_NUM_DOFS + support.iDirection;

		for(size_t i = 0; i < SystemMatrix.Size(); i++) {
			SystemMatrix(i,I) = 0.0;
			SystemMatrix(I,i) = 0.0;
		}

		SystemMatrix(I,I) = 1.0;
		SystemForces(I) = 0.0;
	}

	// Выгрузка перемещений системы в вектор перемещений КЭ
	void PlaceElementGlobalDofs(size_t iElement, ElementVector& elementGlobalDofs)
	{
		Element element = model.Elements(iElement);

		for(size_t i = 0; i < ELEMENT_MATRIX_SIZE; i++) {
			size_t iLocalNode = elementGlobalDofs.LocalNodeId(i);
			size_t iDirection = elementGlobalDofs.DirectionId(i);
			size_t iGlobalNode = element.iNode[iLocalNode];

			elementGlobalDofs(i) = SystemDofs(iGlobalNode, iDirection);
		}
	}
};

// Файл исходных данных
class TextFile {
public:
	TextFile(const char* filename):
		file(filename)
	{
		ReadNextLine();
	}

	// Построчный анализ.
	// Если строка опознана - читается следующая строка и возвращается true.
	// Если строка не опознана - возвращается false и строка остается старой для анализа следующей функцией.

	// Чтение строки заголовка секции
	bool ReadHeaderLine(const char* name)
	{
		if(line != name)
			return false;
		ReadNextLine();
		return true;
	}

	// Чтение строки комментария секции
	bool ReadCommentLine()
	{
		if(line.front() != '['  ||  line.back() != ']')
			return false;
		ReadNextLine();
		return true;
	}

	// Чтение строки узла
	bool ReadNodeLine(Node& node)
	{
		stringstream sLine(line);
		if(!(sLine >> node.x))
			return false;
		if(!(sLine >> node.y))
			return false;
		if(!sLine.eof())
			return false;
		ReadNextLine();
		return true;
	}

	// Чтение строки материала
	bool ReadMaterialLine(Material& material)
	{
		stringstream sLine(line);
		if(!(sLine >> material.A))
			return false;
		if(!(sLine >> material.W))
			return false;
		if(!(sLine >> material.I))
			return false;
		if(!(sLine >> material.E))
			return false;
		if(!sLine.eof())
			return false;
		ReadNextLine();
		return true;
	}

	// Чтение строки элемента
	bool ReadElementLine(Element& element)
	{
		stringstream sLine(line);
		if(!(sLine >> element.iMaterial))
			return false;
		if(!(sLine >> element.iNode[0]))
			return false;
		if(!(sLine >> element.iNode[1]))
			return false;
		if(!sLine.eof())
			return false;
		ReadNextLine();
		return true;
	}

	// Чтение строки связи
	bool ReadSupportLine(Support& support)
	{
		stringstream sLine(line);
		if(!(sLine >> support.iNode))
			return false;
		if(!(sLine >> support.iDirection))
			return false;
		if(!sLine.eof())
			return false;
		ReadNextLine();
		return true;
	}

	// Чтение строки нагрузки
	bool ReadForceLine(Force& force)
	{
		stringstream sLine(line);
		if(!(sLine >> force.iNode))
			return false;
		if(!(sLine >> force.iDirection))
			return false;
		if(!(sLine >> force.value))
			return false;
		if(!sLine.eof())
			return false;

		ReadNextLine();
		return true;
	}

	// Чтение пустой строки
	bool ReadEmptyLine()
	{
		if(!line.empty())
			return false;
		ReadNextLine();
		return true;
	}

	// Контроль достижения конца файла
	bool ReadEof()
	{
		return file.eof();
	}
private:
	ifstream file;
	string line;
private:
	void ReadNextLine()
	{
		line.clear();
		getline(file, line);
	}
};

// Класс чтения файла исходных данных
class ModelImporter {
public:
	ModelImporter(TextFile& textFile, Model& model):
		textFile(textFile), model(model), lineNumber(1)
	{

	}

	bool Import()
	{
		return	ImportNodes() &&
			ImportMaterials() &&
			ImportElements() &&
			ImportSupports() &&
			ImportForces() &&
			End();
	}

	int ShowError()
	{
		cout << "Error in line " << lineNumber << endl;
		return 0;
	}
private:
	TextFile& textFile;
	Model& model;
	size_t lineNumber;
private:
	bool ImportNodes()
	{
		// Чтение заголовка
		if(!textFile.ReadHeaderLine("NODES"))
			return false;
		lineNumber++;

		// Чтение комментария
		if(!textFile.ReadCommentLine())
			return false;
		lineNumber++;

		// Временный вектор для накопления прочитанных узлов
		list<Node> nodes;

		// Чтение узлов
		Node node;
		while(textFile.ReadNodeLine(node)) {
			nodes.push_back(node);
			lineNumber++;
		}

		// Контроль наличия пустой строки в конце секции
		if(!textFile.ReadEmptyLine())
			return false;
		lineNumber++;

		// Если узлов нет - ошибка
		if(nodes.size() == 0)
			return false;

		// Загрузка прочитанных узлов в модель
		model.Nodes.Allocate(nodes.size());
		size_t i = 0;
		for(auto node = nodes.begin(); node != nodes.end(); node++)
			model.Nodes(i++) = *node;

		return true;
	}

	bool ImportMaterials()
	{
		// Чтение заголовка
		if(!textFile.ReadHeaderLine("MATERIALS"))
			return false;
		lineNumber++;

		// Чтение комментария
		if(!textFile.ReadCommentLine())
			return false;
		lineNumber++;

		// Временный вектор для накопления прочитанных материалов
		list<Material> materials;

		// Чтение материалов
		Material material;
		while(textFile.ReadMaterialLine(material)) {
			materials.push_back(material);
			lineNumber++;
		}

		// Контроль наличия пустой строки в конце секции
		if(!textFile.ReadEmptyLine())
			return false;
		lineNumber++;

		// Если материалов нет - ошибка
		if(materials.size() == 0)
			return false;

		// Загрузка прочитанных материалов в модель
		model.Materials.Allocate(materials.size());
		size_t i = 0;
		for(auto material = materials.begin(); material != materials.end(); material++)
			model.Materials(i++) = *material;

		return true;
	}

	bool ImportElements()
	{
		// Чтение заголовка
		if(!textFile.ReadHeaderLine("ELEMENTS"))
			return false;
		lineNumber++;

		// Чтение комментария
		if(!textFile.ReadCommentLine())
			return false;
		lineNumber++;

		// Временный вектор для накопления прочитанных элементов
		list<Element> elements;

		// Чтение и контроль элементов
		Element element;
		while(textFile.ReadElementLine(element)) {
			if(element.iMaterial >= model.Materials.Size())
				return false;
			if(element.iNode[0] >= model.Nodes.Size())
				return false;
			if(element.iNode[1] >= model.Nodes.Size())
				return false;

			elements.push_back(element);
			lineNumber++;
		}

		// Контроль наличия пустой строки в конце секции
		if(!textFile.ReadEmptyLine())
			return false;
		lineNumber++;

		// Если элементов нет - ошибка
		if(elements.size() == 0)
			return false;

		// Загрузка прочитанных элементов в модель
		model.Elements.Allocate(elements.size());
		size_t i = 0;
		for(auto element = elements.begin(); element != elements.end(); element++)
			model.Elements(i++) = *element;

		return true;
	}

	bool ImportSupports()
	{
		// Чтение заголовка
		if(!textFile.ReadHeaderLine("SUPPORTS"))
			return false;
		lineNumber++;

		// Чтение комментария
		if(!textFile.ReadCommentLine())
			return false;
		lineNumber++;

		// Временный вектор для накопления прочитанных связей
		list<Support> supports;

		// Чтение и контроль связей
		Support support;
		while(textFile.ReadSupportLine(support)) {
			if(support.iDirection >= NODE_NUM_DOFS)
				return false;
			if(support.iNode >= model.Nodes.Size())
				return false;

			supports.push_back(support);
			lineNumber++;
		}

		// Контроль наличия пустой строки в конце секции
		if(!textFile.ReadEmptyLine())
			return false;
		lineNumber++;

		// Если связей нет - ошибка
		if(supports.size() == 0)
			return false;

		// Загрузка прочитанных связей в модель
		model.Supports.Allocate(supports.size());
		size_t i = 0;
		for(auto support = supports.begin(); support != supports.end(); support++)
			model.Supports(i++) = *support;

		return true;
	}

	bool ImportForces()
	{
		// Чтение заголовка
		if(!textFile.ReadHeaderLine("FORCES"))
			return false;
		lineNumber++;

		// Чтение комментария
		if(!textFile.ReadCommentLine())
			return false;
		lineNumber++;

		// Временный вектор для накопления прочитанных нагрузок
		list<Force> forces;

		// Чтение и контроль нагрузок
		Force force;
		while(textFile.ReadForceLine(force)) {
			if(force.iNode >= model.Nodes.Size())
				return false;
			if(force.iDirection >= NODE_NUM_DOFS)
				return false;

			forces.push_back(force);
			lineNumber++;
		}

		// Контроль наличия пустой строки в конце секции
//		if(!textFile.ReadEmptyLine())
//			return false;
		lineNumber++;

		// Если нагрузок нет - ошибка
		if(forces.size() == 0)
			return false;

		// Загрузка прочитанных нагрузок в модель
		model.Forces.Allocate(forces.size());
		size_t i = 0;
		for(auto force = forces.begin(); force != forces.end(); force++)
			model.Forces(i++) = *force;

		return true;
	}

	bool End()
	{
		return textFile.ReadEof();
	}
};

// Файл результатов расчета для GMSH
class PosFile {
public:
	PosFile(const char* filename): file(filename) {
		file << "ViewIndex = -1;" << endl;
	}

	// Запись заголовка секции вида
	void WriteViewStart(const char* name)
	{
		file << "ViewIndex++;" << endl;
		file << "View \"" << name << "\"{" << endl;
	}

	// Запись точки числового поля
	void WriteScalarPoint(double x, double y, double value)
	{
		double z = 0.0;
		file	<< "\tSP ("
			<< x << ","
			<< y << ","
			<< z << ") {"
			<< value << "};" << endl;
	}

	// Запись линии числового поля с постоянным значением
	void WriteScalarLine(double x1, double y1, double x2, double y2, double value)
	{
		double z1 = 0.0;
		double z2 = 0.0;
		file	<< "\tSL ("
			<< x1 << ","
			<< y1 << ","
			<< z1 << ","
			<< x2 << ","
			<< y2 << ","
			<< z2 << ") {"
			<< value << ","
			<< value << "};" << endl;
	}

	// Запись линии числового поля с переменным значением
	void WriteScalarLine(double x1, double y1, double x2, double y2, double value1, double value2)
	{
		double z1 = 0.0;
		double z2 = 0.0;
		file	<< "\tSL ("
			<< x1 << ","
			<< y1 << ","
			<< z1 << ","
			<< x2 << ","
			<< y2 << ","
			<< z2 << ") {"
			<< value1 << ","
			<< value2 << "};" << endl;
	}

	// Запись точки векторного поля
	void WriteVectorPoint(double x, double y, double vx, double vy, double vz)
	{
		double z = 0.0;
		file	<< "\tVP ("
			<< x << ","
			<< y << ","
			<< z << ") {"
			<< vx << ","
			<< vy << ","
			<< vz << "};" << endl;
	}

	// Запись завершение секции вида
	void WriteViewEnd()
	{
		file << "};" << endl;
	}

	// Запись настроек вида
	void WriteViewOption_IntervalsType(int intervalsType)
	{
		file << "View[ViewIndex].IntervalsType = " << intervalsType << ";" << endl;
	}

	void WriteViewOption_NbIso(int nbIso)
	{
		file << "View[ViewIndex].NbIso = " << nbIso << ";" << endl;
	}

	void WriteViewOption_LineWidth(double lineWidth)
	{
		file << "View[ViewIndex].LineWidth = " << lineWidth << ";" << endl;
	}

	void WriteViewOption_Format(const char* format)
	{
		file << "View[ViewIndex].Format = \"" << format << "\";" << endl;
	}

	void WriteViewOption_Visible(bool visible)
	{
		file << "View[ViewIndex].Visible = " << visible << ";" << endl;
	}

	void WriteViewOption_ColormapNumber(bool colormapNumber)
	{
		file << "View[ViewIndex].ColormapNumber = " << colormapNumber << ";" << endl;
	}

	void WriteViewOption_ShowScale(bool showScale)
	{
		file << "View[ViewIndex].ShowScale = " << showScale << ";" << endl;
	}

	void WriteViewOption_VectorType(int vectorType)
	{
		file << "View[ViewIndex].VectorType = " << vectorType << ";" << endl;
	}

	void WriteViewOption_CenterGlyphs(int centerGlyphs)
	{
		file << "View[ViewIndex].CenterGlyphs = " << centerGlyphs << ";" << endl;
	}
private:
	ofstream file;
};

// Класс экспорта данных модели в файл GMSH
class ModelExporter {
public:
	ModelExporter(PosFile& posFile, Model& model):
		posFile(posFile), model(model)
	{

	}

	void Export()
	{
		ExportGeometry();
		ExportNodesNumbers();
		ExportElementNumbers();
		ExportMaterials();
		ExportForces();
		ExportSupports();
	}
private:
	PosFile& posFile;
	Model& model;
private:
	void ExportGeometry()
	{
		posFile.WriteViewStart("Geometry");

		for(size_t i = 0; i < model.Elements.Size(); i++) {
			Element element = model.Elements(i);
			Node node0 = model.Nodes(element.iNode[0]);
			Node node1 = model.Nodes(element.iNode[1]);
			posFile.WriteScalarLine(node0.x, node0.y, node1.x, node1.y, i);
		}

		posFile.WriteViewEnd();

		posFile.WriteViewOption_ColormapNumber(0);
		posFile.WriteViewOption_ShowScale(false);
	}

	void ExportNodesNumbers()
	{
		posFile.WriteViewStart("Nodes numbers");

		for(size_t i = 0; i < model.Nodes.Size(); i++) {
			Node node = model.Nodes(i);
			posFile.WriteScalarPoint(node.x, node.y, i);
		}

		posFile.WriteViewEnd();

		posFile.WriteViewOption_ColormapNumber(0);
		posFile.WriteViewOption_IntervalsType(4);
		posFile.WriteViewOption_ShowScale(false);
		posFile.WriteViewOption_Visible(false);
	}

	void ExportElementNumbers()
	{
		posFile.WriteViewStart("Element numbers");

		for(size_t i = 0; i < model.Elements.Size(); i++) {
			Element element = model.Elements(i);
			Node node0 = model.Nodes(element.iNode[0]);
			Node node1 = model.Nodes(element.iNode[1]);
			posFile.WriteScalarLine(node0.x, node0.y, node1.x, node1.y, i);
		}

		posFile.WriteViewEnd();

		posFile.WriteViewOption_ColormapNumber(0);
		posFile.WriteViewOption_IntervalsType(4);
		posFile.WriteViewOption_ShowScale(false);
		posFile.WriteViewOption_Visible(false);
	}

	void ExportMaterials()
	{
		posFile.WriteViewStart("Materials");

		for(size_t i = 0; i < model.Elements.Size(); i++) {
			Element element = model.Elements(i);
			Node node0 = model.Nodes(element.iNode[0]);
			Node node1 = model.Nodes(element.iNode[1]);
			posFile.WriteScalarLine(node0.x, node0.y, node1.x, node1.y, element.iMaterial);
		}

		posFile.WriteViewEnd();

		posFile.WriteViewOption_IntervalsType(3);
		posFile.WriteViewOption_NbIso(model.Materials.Size());
		posFile.WriteViewOption_LineWidth(10);
		posFile.WriteViewOption_Format("%.0f");
		posFile.WriteViewOption_Visible(false);
	}

	void ExportForces()
	{
		posFile.WriteViewStart("Forces");

		for(size_t i = 0; i < model.Forces.Size(); i++) {
			Force force = model.Forces(i);
			Node node = model.Nodes(force.iNode);

			double v[] = {0,0,0};
			v[force.iDirection] = force.value;

			posFile.WriteVectorPoint(node.x, node.y, v[0], v[1], v[2]);
		}

		posFile.WriteViewEnd();
		posFile.WriteViewOption_ColormapNumber(1);
		posFile.WriteViewOption_ShowScale(false);
	}

	void ExportSupports()
	{
		posFile.WriteViewStart("Supports");

		for(size_t i = 0; i < model.Supports.Size(); i++) {
			Support support = model.Supports(i);
			Node node = model.Nodes(support.iNode);

			double v[] = {0,0,0};
			v[support.iDirection] = 1.0;

			posFile.WriteVectorPoint(node.x, node.y, v[0], v[1], v[2]);
		}

		posFile.WriteViewEnd();
		posFile.WriteViewOption_LineWidth(10);
		posFile.WriteViewOption_ColormapNumber(0);
		posFile.WriteViewOption_VectorType(1);
		posFile.WriteViewOption_CenterGlyphs(1);
		posFile.WriteViewOption_ShowScale(false);
	}
};

// Класс экспорта данных элементов в файл GMSH
class ElementsDataExporter {
public:
	ElementsDataExporter(PosFile& posFile, Model& model, ElementsData& elementsData):
		posFile(posFile), model(model), elementsData(elementsData)
	{

	}

	void Export()
	{
		ExportElementResults("N",     elementsData.LongitudinalForces);
		ExportElementResults("Q",     elementsData.ShearForces);
		ExportElementResults("M",     elementsData.BendingMoments);
		ExportElementResults("Sigma", elementsData.Stresses);
	}
private:
	PosFile& posFile;
	Model& model;
	ElementsData& elementsData;
private:
	// Запись вектора результатов в расчетных сечениях КЭ
	void ExportElementResults(const char* name, Array<ElementResult>& elementResults)
	{
		posFile.WriteViewStart(name);

		for(size_t i = 0; i < model.Elements.Size(); i++) {
			ElementResult value = elementResults(i);

			Element element = model.Elements(i);
			Node node0 = model.Nodes(element.iNode[0]);
			Node node1 = model.Nodes(element.iNode[1]);

			// Координаты промежуточных сечений КЭ
			double x[ELEMENT_NUM_SECTIONS];
			double y[ELEMENT_NUM_SECTIONS];

			// Вычисление координат промежуточных сечений
			const int N = ELEMENT_NUM_SECTIONS - 1;
			for(size_t n = 0; n < ELEMENT_NUM_SECTIONS; n++) {
				double k0 = double(N - n) / double(N);
				double k1 = double(n)     / double(N);

				x[n] = k0 * node0.x + k1 * node1.x;
				y[n] = k0 * node0.y + k1 * node1.y;
			}

			// Запись значений по всем промежуточным сечениям
			for(size_t n = 0; n < N; n++) {
				posFile.WriteScalarLine(x[n], y[n], x[n+1], y[n+1], value(n), value(n+1));
			}
		}

		posFile.WriteViewEnd();

		posFile.WriteViewOption_ColormapNumber(1);
		posFile.WriteViewOption_LineWidth(10);
		posFile.WriteViewOption_Visible(false);
	}
};

// Отладочные функции вывода на экран матриц и векторов
ostream& operator<< (ostream& cout, ElementMatrix& matrix)
{
	cout << "ElementMatrix:" << endl;
	for(size_t i = 0; i < ELEMENT_MATRIX_SIZE; i++) {
		for(size_t j = 0; j < ELEMENT_MATRIX_SIZE; j++)
			cout << " " << matrix(i,j);
		cout << endl;
	}
	return cout << endl;
}

ostream& operator<< (ostream& cout, TransformationMatrix& matrix)
{
	cout << "TransformationMatrix:" << endl;
	for(size_t i = 0; i < ELEMENT_MATRIX_SIZE; i++) {
		for(size_t j = 0; j < ELEMENT_MATRIX_SIZE; j++)
			cout << " " << matrix(i,j);
		cout << endl;
	}
	return cout << endl;
}

ostream& operator<< (ostream& cout, ElementVector& vector)
{
	cout << "ElementVector:" << endl;
	for(size_t i = 0; i < ELEMENT_MATRIX_SIZE; i++)
		cout << " " << vector(i) << endl;
	return cout << endl;
}

ostream& operator<< (ostream& cout, ElementResult& vector)
{
	cout << "ElementResult:" << endl;
	for(size_t i = 0; i < ELEMENT_NUM_SECTIONS; i++)
		cout << " " << vector(i) << endl;
	return cout << endl;
}

ostream& operator<< (ostream& cout, SystemMatrix& matrix)
{
	cout << "SystemMatrix:" << endl;
	for(size_t i = 0; i < matrix.Size(); i++) {
		for(size_t j = 0; j < matrix.Size(); j++)
			cout << " " << matrix(i,j);
		cout << endl;
	}
	return cout << endl;
}

ostream& operator<< (ostream& cout, SystemVector& vector)
{
	cout << "SystemVector:" << endl;
	for(size_t i = 0; i < vector.Size(); i++)
		cout << " " << vector(i) << endl;
	return cout << endl;
}

int main(int argc, char* argv[])
{
	Model model;

	// Чтение исходных данных
	TextFile textFile("FE.txt");
	ModelImporter modelImporter(textFile, model);
	if(!modelImporter.Import())
		return modelImporter.ShowError();

	// Расчет
	ElementsData elementsData(model);

	elementsData.CalcLocalMatrices();
	elementsData.CalcTransformationMatrices();
	elementsData.CalcGlobalMatrices();

	Solver solver(model, elementsData);

	solver.PlaceElementGlobalMatrices();
	solver.PlaceForces();
	solver.PlaceSupports();
	solver.SolveGauss();
	solver.PlaceElementGlobalDofs();

	elementsData.CalcLocalDofs();
	elementsData.CalcLocalReactions();
	elementsData.CalcLongitudinalForces();
	elementsData.CalcShearForces();
	elementsData.CalcBendingMoments();
	elementsData.CalcStresses();

	// Запись результатов
	PosFile posFile("FE.pos");
	ModelExporter modelExporter(posFile, model);
	modelExporter.Export();

	ElementsDataExporter elementsDataExporter(posFile, model, elementsData);
	elementsDataExporter.Export();

	/* Отладочный вывод матриц на экран
	cout << "SystemMatrix" << endl;
	cout << solver.SystemMatrix;

	cout << "SystemDofs" << endl;
	cout << solver.SystemDofs;

	cout << "GlobalDofs" << endl;
	for(size_t i = 0; i < elementsData.GlobalDofs.Size(); i++)
		cout << elementsData.GlobalDofs(i);

	cout << "LocalDofs" << endl;
	for(size_t i = 0; i < elementsData.LocalDofs.Size(); i++)
		cout << elementsData.LocalDofs(i);

	cout << "LocalReactions" << endl;
	for(size_t i = 0; i < elementsData.LocalReactions.Size(); i++)
		cout << elementsData.LocalReactions(i);

	cout << "LongitudinalForces" << endl;
	for(size_t i = 0; i < elementsData.LongitudinalForces.Size(); i++)
		cout << elementsData.LongitudinalForces(i);

	cout << "ShearForces" << endl;
	for(size_t i = 0; i < elementsData.ShearForces.Size(); i++)
		cout << elementsData.ShearForces(i);

	cout << "BendingMoments" << endl;
	for(size_t i = 0; i < elementsData.BendingMoments.Size(); i++)
		cout << elementsData.BendingMoments(i);

	cout << "Stresses" << endl;
	for(size_t i = 0; i < elementsData.Stresses.Size(); i++)
		cout << elementsData.Stresses(i);
	*/
}

Программа читает файл входных данных текстового формата (нумерация узлов, материалов, элементов - от нуля):

Код:

[Выделить все]

NODES
[X Y]
0 0
0 1
1 1

MATERIALS
[A W I E]
1 1 1 1
1 1 1 1

ELEMENTS
[Material Node1 Node2]
0 0 1
1 1 2

SUPPORTS
[Node Direction]
0 0
0 1
0 2

FORCES
[Node Direction Value]
2 0 0.5
2 1 1.0

Записывает результаты в файл FE.pos для просмотра в GMSH

Все дубовое - плотные матрицы, неучет симметрии и т.п., в расчете на задачи в десяток-другой элементов. Пока сделано то, что есть в каждой книге по МКЭ - решение статической задачи. Теперь надо прикрутить жесткие вставки для задания начальных искривлений, вывести матрицы жесткости реакций и дополнительных реакций, добавить геомнелин и поиск минимума напряжений по вектору искривлений.

Во вложении - скомпилированная в VS версия с тестовым файлом.

[румата]

Сообщение от Нубий-IV Заготовка для обработки напильником.

МКЭ решатель в 1500 строк кода?

[Бахил]

Сообщение от румата МКЭ решатель в 1500 строк кода?

Действительно. Понаплодил ненужных структур и классов. Вполне можно было в 100 строк уложитбся.
Ну и где там геометрическая нелинейность?

[Нубий-IV]

Для упрощения счета все параметры (силы, длины, жесткости) единичные.

• Статический расчет
  [K]z=F
  
  
  
  
  
  
  Прогиб 0.3333
  Момент в заделке M=1
• Расчет на устойчивость
  [K+NG]z = 0
  
  
  
  
  Ncr=2.486
• Геоменлин
  [K+NG]z = F + [K]z0
  
  
  
  
  
  
  
  
  
  
  
  Прогиб 1.1127
  Момент в заделке M=2.1127
• Ручной счет
  
  
  
  Прогиб 1.116
  Момент в заделке M=2.116

У этих формул интересные последствия:

• При фиксированной N начальное искривление - всего лишь вид нагрузки.
• Следствие - напряжение в элементе зависит от начального искривления линейно.
• Чтобы найти наихудшую форму начального искривления, достаточно для каждого узла с возможным иск взять знак, увеличивающий напряжение. Эта форма даст наибольшее напряжение при заданной нагрузке.
• Работает принцип "максимальной жадности" - величина искривления в узле максимальна в пределах допуска.
• Расчет аналогичен загружению линии влияния, только матрица системы равна не [K], а [K-NG].
• Такой расчет можно использовать для контроля, когда нагрузка уже найдена из какого-то другого расчета: можно создать схему с опасным для нее начальным искривлением и проверить в честном нелинейном счете.
• При других нагрузках меняется матрица системы [K-NG], и искривление может стать другим. Этот эффект виден в чудо-раме, если смотреть загружения по шагам. Сначала побеждает внешняя нагрузка, потом - форма потери устойчивости.
• Чтобы решить честную задачу определения N, надо повторять расчет при других значениях нагрузки. Но это уже задача с одним неизвестным N, а не перебор всех возможных искривлений. Комбинаторный взрыв отменяется. Одно нелинейное уравнение - это уже не то, что нелинейный поиск по всем степеням свободы системы.

Сообщение от румата МКЭ решатель в 1500 строк кода

В память об одной теме, которая померла, не родившись .

Сообщение от Бахил можно было в 100 строк уложитбся. Ну и где там геометрическая нелинейность?

Работающая - посте 405. А в тот, кто пишет МКЭ за 100 строк (это же всего час - другой, да?), может по формулам выше свалять свою версию и показать на форуме.

[Бахил]

Ошибки в матрицах.
Если справа жёсткая заделка, то всего 3 неизвестных.
Если справа единичные связи КЖ, то последние 3 строки неверны.
И продольная свобода в данном примере совершенно лишняя.
Offtop: Про 100 строк, это сарказм на пост руматы (1500).

[Нубий-IV]

Сообщение от Бахил всего 3 неизвестных

Это один из способов учета связей - зануление всех элементов, кроме диагонали, где ставятся единицы; вектор нагрузок в строках тоже зануляется. Это гарантирует нулевые ответы в таких строках, что соответствует связям. Плюс не надо удалять индексы и выносить себе мозг вычислением смещений после перенумерации из-за вычеркнутых строк. И вектор неизвестных сразу полный, не надо дописывать обратно вычеркнутые строки, нули там появятся в ходе решения системы, а не от лишних циклов с нулями. В моей программе связи ставятся так же.

Сообщение от Бахил продольная свобода в данном примере совершенно лишняя

Что значит "лишняя"?! Это продольная сила, которая и создает нелинейность. Ради нее все и затевалось. И она дала почти точное сопроматовское решение задачи с множителем 1/(1-N/Ncr). Причем без всяких там шаговых методов, в один Гаусс. Это простейшая задача на устойчивую прочность с учетом начального смещения и боковой нагрузки.

[Yu Mo]

Спустимся на землю.
Десятые лировцы показали в YouTube свою небольшую, но полезную фишку по определению расчётных длин стоек с неполными связями. Никаких революций. Понятные рассуждения, основанные на принятой в большинстве учебников и справочников теории, привели к конкретному и ожидаемому результату.
Думаю, заценить это могут металлисты, работающие в отраслевых институтах и фирмах, там, где первую скрипку играют технологи. А они, в общем-то, враждебно относятся не только к связям между стойками, но и к самим стойкам. Ну, со стойками они как-то уж смирились, а вот насчёт связей, возможны дискуссии. Однако же, бывая на наших объектах, я, увы, вынужден был иногда признать их правоту.

https://www.youtube.com/watch?v=sBTGcSlsm_I

P.S. Если из их сообщения слить всю воду, то в сухом остатке останутся три таблицы которые можно скопировать и пользоваться ими при проектировании. При этом верхнее ограничение коэффициента расчётной длины, упоминаемое в #33, #38, на совести проектировщика.

[Бам]

Сообщение от Yu Mo Десятые лировцы показали в YouTube свою небольшую, но полезную фишку по определению расчётных длин стоек с неполными связями. Никаких революций. Понятные рассуждения, основанные на принятой в большинстве учебников и справочников теории, привели к конкретному и ожидаемому результату.

А экспертиза такие расчеты уже принимает, лировцы не говорят? Серьезная промка это еще и главгосэкспертиза зачастую....

[Бахил]

Сообщение от Нубий-IV Это один из способов учета связей

Весьма нерациональный. Всё занулять необязательно. Достаточно диагональные элементы увеличить на порядок.

Сообщение от Нубий-IV Что значит "лишняя"?!

То и значит - она никак не участвует в определении критической силы.

Сообщение от Нубий-IV Причем без всяких там шаговых методов

Ну так линейная алгебра на то и "линейная", что не требует никаких шаговых методов.

Сообщение от Нубий-IV [K+NG]z = 0