Реклама i
Стань нашим партнером - зарабатывай с нами!
| Правила | Регистрация | Пользователи | Поиск | Сообщения за день | Все разделы прочитаны |  Справка по форуму | Файлообменник |

Вернуться   Форум DWG.RU > Архитектура и Строительство > Конструкции зданий и сооружений > Металлические конструкции > Методы определения расчетных длин пригодных для расчетов на устойчивость по СП 16.13330. Ищем, делимся, обсуждаем.

Методы определения расчетных длин пригодных для расчетов на устойчивость по СП 16.13330. Ищем, делимся, обсуждаем.

Ответ
Поиск в этой теме
Непрочитано 09.09.2022, 22:37 1 |
Методы определения расчетных длин пригодных для расчетов на устойчивость по СП 16.13330. Ищем, делимся, обсуждаем.
румата
 
Регистрация: 06.04.2015
Сообщений: 2,673

Эта тема создана с целью поиска и сбора информации по любым из существующих способов и методам определения расчетных длин, а также с целью их обсуждения, критики или формулировок своих собственный методик применительно к расчетам на устойчивость и проверке по предельной гибкости согласно норм РФ и, возможно, других стран.
Просмотров: 120413
 
Непрочитано 26.10.2022, 16:36
1 | 1 #461
Бам


 
Регистрация: 06.09.2022
Сообщений: 974


Я так понимаю, по всей стране всё расчеты сжатых рам нужно немедленно прервать, пока ваши исследования не будут закончены? Нас ждут сенсационные результаты?
Бам вне форума  
 
Непрочитано 27.10.2022, 13:11
#462
Ильнур

КМ (+КМД), КЖ (КЖФ)
 
Регистрация: 30.05.2007
Далече
Сообщений: 25,086


Цитата:
Сообщение от Бахил Посмотреть сообщение
... А длины то будут?
Это несходящийся процесс.
Цитата:
Нас ждут сенсационные результаты?
Да, будет показано, что правильная расчетная длина - величина неопределяемая в принципе, если стержней более 1 шт.
__________________
Воскресе
Ильнур вне форума  
 
Непрочитано 27.10.2022, 13:28
#463
Бахил

?
 
Регистрация: 17.06.2014
Царицын
Сообщений: 12,203


Цитата:
Сообщение от Ильнур Посмотреть сообщение
правильная расчетная длина - величина неопределяемая в принципе, если стержней более 1 шт.
Да! Это надо отлить в бронзе и поместить на фронтоны всех строительных вузов.
__________________
Не откладывайте на завтра! Положите на всё уже сегодня.(с)
Бахил вне форума  
 
Непрочитано 27.10.2022, 14:39
#464
crossing


 
Регистрация: 19.07.2018
Сообщений: 3,081


Надо вот этого вот румату забанить на всегда. Как забанили этого семсонофф и онже илюха. И румата из тогоже набора.
crossing вне форума  
 
Автор темы   Непрочитано 27.10.2022, 15:18
#465
румата


 
Регистрация: 06.04.2015
Сообщений: 2,673


Цитата:
Сообщение от crossing Посмотреть сообщение
Надо вот этого вот румату забанить на всегда.
Offtop: Интересно, а за что надо банить румату?
румата вне форума  
 
Непрочитано 27.10.2022, 15:24
#466
Бам


 
Регистрация: 06.09.2022
Сообщений: 974


Offtop:
Цитата:
Сообщение от румата Посмотреть сообщение
Offtop: Интересно, а за что надо банить румату?
"А вы не любите пролетариат..."(с)
Бам вне форума  
 
Непрочитано 27.10.2022, 17:58
#467
crossing


 
Регистрация: 19.07.2018
Сообщений: 3,081


Цитата:
Сообщение от румата Посмотреть сообщение
Offtop: Интересно, а за что надо банить румату?
За то, что не понимаешь о чём спрашиваешь. Вот ты тему создал, а у тебя в мозгах контекста нет.
Ты через полгода ещё такуюже сосдаш.
crossing вне форума  
 
Автор темы   Непрочитано 27.10.2022, 18:12
#468
румата


 
Регистрация: 06.04.2015
Сообщений: 2,673


Цитата:
Сообщение от crossing Посмотреть сообщение
За то, что не понимаешь о чём спрашиваешь.
В теме вообще нет вопроса. Эта тема для сбора информации и поиска способов и методов определния расчетных длин пригодных для расчетов по СП16.13330.
Цитата:
Сообщение от crossing Посмотреть сообщение
Вот ты тему создал, а у тебя в мозгах контекста нет.
Снова интересно. Что ты подразумеваешь под контекстом, которого у меня нет в мозгах?
Цитата:
Сообщение от crossing Посмотреть сообщение
Ты через полгода ещё такуюже сосдаш.
Зачем? Мне этой вполне достатотчно на всю оставшуюся жизнь.
румата вне форума  
 
Непрочитано 01.11.2022, 05:29
#469
Нубий-IV

Инженер-философ
 
Регистрация: 24.04.2019
Хабаровск
Сообщений: 1,867


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

Например, для стержня Эйлера при учете пары первых форм:
  • В первом узле коэффициенты формы связаны условием вписаться в допуск:

  • Аналогично, во втором узле:

  • Если перебрать все узлы - получится область допустимых значений.
  • Если на эту область нанести изолинии постоянных N, можно найти наихудшую комбинацию:
    Сейчас изолинии почти вертикальные, но это особенность одиночного стержня, где вторая критическая в разы выше первой. В рамах, если соседние формы почти равны, линии наклонятся ближе к 45°; а если и соотношение форм удачно ляжет - то наклонятся еще сильнее. Линии не параллельны, они постепенно выходят на вертикаль при приближении к Ncr1, или к фиксированному углу, если RA<Ncr1.
  • Сразу можно отбраковать только случай двух отрицательных коэффициентов (он означает максимальное напряжение в противоположной точке сечения).
    В первой четверти решение есть всегда. В четвертях с разными знаками возможны области, где решения нет. Но в целом при других соотношениях форм и критических сил возможны и варианты, когда надо выбирать единственную форму, и варианты, когда хуже оказываются комбинации. Например, чудо-рама - случай выбора второй формы, а Эйлер с постоянным допуском - случай выбора комбинации.
  • С учетом того, что зависимость N от коэффициентов формы нелинейная, получается задача нелинейного программирования. Сходу можно только сказать, что минимум попадает в узел области, потому что изолинии N прямые, то есть достаточно перебора конечного числа точек.
  • Когда коэффициенты форм выбраны, поиск N - это поиск наименьшего корня полинома. Каждая форма добавляет единицу к степени уравнения. При одной форме уравнение квадратное, при двух - кубическое и т.д.
  • Кто может сказать - проще это, чем поиск минимума среди всех возможных искривлений через последовательное решение системы уравнений MKЭ?
  • В случае с Эйлером подсказку дает перемножение эпюр, это намного проще, чем искать минимум бесконечного числа форм. Для рамы возможны какие-нибудь аналогичные соображения?
Цитата:
Сообщение от Бам Посмотреть сообщение
по всей стране всё расчеты сжатых рам нужно немедленно прервать, пока ваши исследования не будут закончены?
Риторические вопросы слишком сложны, их надо постить в раздел форума "философия". А продемонстрировать острый ум можно намного проще. Например, предоставить готовый ответ на вопрос темы. Или хотя бы ткнуть в какую-нибудь ошибку (в посте 457 есть минимум одна; не буду удалять, пусть лежит).
Цитата:
Сообщение от Ильнур Посмотреть сообщение
расчетная длина - величина неопределяемая в принципе, если стержней более 1 шт
Цитата:
Сообщение от Бахил Посмотреть сообщение
Да! Это надо отлить в бронзе и поместить на фронтоны всех строительных вузов.
Цитата:
Сообщение от румата Посмотреть сообщение
crossing Ты через полгода ещё такуюже сосдаш.
румата Зачем?
Если у вас нету темы,
В нее не нагадит сосед.
И не обсуждать вам проблемы
С тем, для кого,
С тем, для кого,
С тем, для кого их нет.
Миниатюры
Нажмите на изображение для увеличения
Название: 01.png
Просмотров: 662
Размер:	80.9 Кб
ID:	250912  Нажмите на изображение для увеличения
Название: 02.png
Просмотров: 673
Размер:	79.9 Кб
ID:	250913  Нажмите на изображение для увеличения
Название: 03.png
Просмотров: 671
Размер:	44.4 Кб
ID:	250914  Нажмите на изображение для увеличения
Название: 04.jpg
Просмотров: 667
Размер:	156.6 Кб
ID:	250915  Нажмите на изображение для увеличения
Название: 05.PNG
Просмотров: 646
Размер:	2.1 Кб
ID:	250916  

Нубий-IV вне форума  
 
Непрочитано 01.11.2022, 09:21
#470
Бахил

?
 
Регистрация: 17.06.2014
Царицын
Сообщений: 12,203


"Формы" линейно независимы и ортогональны. Фундаментальное свойство собственных чисел и векторов.
__________________
Не откладывайте на завтра! Положите на всё уже сегодня.(с)
Бахил вне форума  
 
Непрочитано 01.11.2022, 16:34
#471
Нубий-IV

Инженер-философ
 
Регистрация: 24.04.2019
Хабаровск
Сообщений: 1,867


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

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

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

Пока что вижу только решение "в лоб" - свалять свой микро-мкэ для плоских рам, и в нем поиск минимума по начальным искривлениям организовать для схем из СП. Но запрограммить это небыстро будет.
Нубий-IV вне форума  
 
Непрочитано 18.11.2022, 14:41
1 | #472
Нубий-IV

Инженер-философ
 
Регистрация: 24.04.2019
Хабаровск
Сообщений: 1,867


Заготовка для обработки напильником. Исходник 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 версия с тестовым файлом.
Миниатюры
Нажмите на изображение для увеличения
Название: FE.png
Просмотров: 573
Размер:	54.0 Кб
ID:	251335  
Вложения
Тип файла: zip FE.ZIP (23.9 Кб, 16 просмотров)
Нубий-IV вне форума  
 
Автор темы   Непрочитано 21.11.2022, 10:56
#473
румата


 
Регистрация: 06.04.2015
Сообщений: 2,673


Цитата:
Сообщение от Нубий-IV Посмотреть сообщение
Заготовка для обработки напильником.
МКЭ решатель в 1500 строк кода?
румата вне форума  
 
Непрочитано 21.11.2022, 15:20
#474
Бахил

?
 
Регистрация: 17.06.2014
Царицын
Сообщений: 12,203


Цитата:
Сообщение от румата Посмотреть сообщение
МКЭ решатель в 1500 строк кода?
Действительно. Понаплодил ненужных структур и классов. Вполне можно было в 100 строк уложитбся.
Ну и где там геометрическая нелинейность?
__________________
Не откладывайте на завтра! Положите на всё уже сегодня.(с)
Бахил вне форума  
 
Непрочитано 22.11.2022, 08:57
1 | #475
Нубий-IV

Инженер-философ
 
Регистрация: 24.04.2019
Хабаровск
Сообщений: 1,867


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


  • Статический расчет
    [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 строк (это же всего час - другой, да?), может по формулам выше свалять свою версию и показать на форуме.
Миниатюры
Нажмите на изображение для увеличения
Название: Схема.png
Просмотров: 492
Размер:	6.8 Кб
ID:	251388  Нажмите на изображение для увеличения
Название: K.PNG
Просмотров: 499
Размер:	6.3 Кб
ID:	251389  Нажмите на изображение для увеличения
Название: G.PNG
Просмотров: 499
Размер:	5.3 Кб
ID:	251390  
Нубий-IV вне форума  
 
Непрочитано 22.11.2022, 10:24
#476
Бахил

?
 
Регистрация: 17.06.2014
Царицын
Сообщений: 12,203


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

Последний раз редактировалось Бахил, 22.11.2022 в 11:49.
Бахил вне форума  
 
Непрочитано 22.11.2022, 13:21
#477
Нубий-IV

Инженер-философ
 
Регистрация: 24.04.2019
Хабаровск
Сообщений: 1,867


Цитата:
Сообщение от Бахил Посмотреть сообщение
всего 3 неизвестных
Это один из способов учета связей - зануление всех элементов, кроме диагонали, где ставятся единицы; вектор нагрузок в строках тоже зануляется. Это гарантирует нулевые ответы в таких строках, что соответствует связям. Плюс не надо удалять индексы и выносить себе мозг вычислением смещений после перенумерации из-за вычеркнутых строк. И вектор неизвестных сразу полный, не надо дописывать обратно вычеркнутые строки, нули там появятся в ходе решения системы, а не от лишних циклов с нулями. В моей программе связи ставятся так же.
Цитата:
Сообщение от Бахил Посмотреть сообщение
продольная свобода в данном примере совершенно лишняя
Что значит "лишняя"?! Это продольная сила, которая и создает нелинейность. Ради нее все и затевалось. И она дала почти точное сопроматовское решение задачи с множителем 1/(1-N/Ncr). Причем без всяких там шаговых методов, в один Гаусс. Это простейшая задача на устойчивую прочность с учетом начального смещения и боковой нагрузки.
Нубий-IV вне форума  
 
Непрочитано 22.11.2022, 14:45
#478
Yu Mo

Расчёты и проектирование строительных конструкций
 
Регистрация: 04.03.2008
Донецк
Сообщений: 409


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

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

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


 
Регистрация: 06.09.2022
Сообщений: 974


Цитата:
Сообщение от Yu Mo Посмотреть сообщение
Десятые лировцы показали в YouTube свою небольшую, но полезную фишку по определению расчётных длин стоек с неполными связями. Никаких революций. Понятные рассуждения, основанные на принятой в большинстве учебников и справочников теории, привели к конкретному и ожидаемому результату.
А экспертиза такие расчеты уже принимает, лировцы не говорят? Серьезная промка это еще и главгосэкспертиза зачастую....
Бам вне форума  
 
Непрочитано 22.11.2022, 15:17
#480
Бахил

?
 
Регистрация: 17.06.2014
Царицын
Сообщений: 12,203


Цитата:
Сообщение от Нубий-IV Посмотреть сообщение
Это один из способов учета связей
Весьма нерациональный. Всё занулять необязательно. Достаточно диагональные элементы увеличить на порядок.
Цитата:
Сообщение от Нубий-IV Посмотреть сообщение
Что значит "лишняя"?!
То и значит - она никак не участвует в определении критической силы.
Цитата:
Сообщение от Нубий-IV Посмотреть сообщение
Причем без всяких там шаговых методов
Ну так линейная алгебра на то и "линейная", что не требует никаких шаговых методов.

Цитата:
Сообщение от Нубий-IV Посмотреть сообщение
[K+NG]z = 0
__________________
Не откладывайте на завтра! Положите на всё уже сегодня.(с)

Последний раз редактировалось Бахил, 22.11.2022 в 15:23.
Бахил вне форума  
Ответ
Вернуться   Форум DWG.RU > Архитектура и Строительство > Конструкции зданий и сооружений > Металлические конструкции > Методы определения расчетных длин пригодных для расчетов на устойчивость по СП 16.13330. Ищем, делимся, обсуждаем.

Размещение рекламы
Опции темы Поиск в этой теме
Поиск в этой теме:

Расширенный поиск


Похожие темы
Тема Автор Раздел Ответов Последнее сообщение
В поиске сравнительные таблицы старых (СНиП и пр.) и новых нормативных документов (актуализированные СП) Armin Поиск литературы, чертежей, моделей и прочих материалов 19 25.11.2016 08:27
Как трактовать указания СП 16.13330 "Стальные конструкции"? gdenisn Металлические конструкции 41 20.10.2016 06:37
Обязательные и доброволные нормы Aragorn Прочее. Архитектура и строительство 24 15.12.2014 14:08
Расчет ангара в Scad. Вопрос по коэффициентам расчетных длин для связей. TOWER SCAD 9 15.07.2009 07:46
Коэффициенты расчетных длин в постпроцессоре SCAD Pilot729 SCAD 4 25.12.2006 12:36