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

Вернуться   Форум DWG.RU > Программное обеспечение > Программирование > Как создать расчетное программное обеспечение с открытым исходным кодом (конструктивные решения)

Как создать расчетное программное обеспечение с открытым исходным кодом (конструктивные решения)

Ответ
Поиск в этой теме
Непрочитано 24.09.2021, 14:52
Как создать расчетное программное обеспечение с открытым исходным кодом (конструктивные решения)
nickname2019
 
Регистрация: 18.11.2019
Сообщений: 1,037

На мой взгляд, основными проблемами российского рынка расчетного программного обеспечения являются:
- отсутствие нормальной возможности программной автоматизации по решению расчетных задач
(в расчетных программах отсутствует возможность для нормального программирования, т.е. невозможно написать программу для полностью автоматического создания расчетной схемы (нескольких расчетных схем), автоматического выполнения расчета, автоматического получения результатов и их автоматического анализа);
- закрытый исходный код по подбору расчетных параметров несущих элементов
(различные программы дают различные результаты при решении одинаковых задач, сравнение алгоритмов подбора различных между собой невозможно, так как код закрыт, общепринятых и одобренных алгоритмов нет, каждый пользуется своим "черным ящиком", который иногда может выдать ошибочное решение);
- для людей, которые занимаются автоматизацией на Лисп, C# и т.д. отсутствуют инструменты, которые позволяли бы программно "по-простому" вызвать готовую библиотечную функцию (например, по подбору сечения какой-то простой балки непосредственно из графического редактора), что вызывает необходимость вызова отдельной расчетной программы, что серьезно тормозит работу;
- "корявый" интерфейс, ужасно неудобная и медленная работа в существующих российских (и украинских) расчетных программах;
(фактически при наличии нормального графического редактора (autocad, nanocad и т.д.) приходится экспортировать данные в в "корявый" редактор расчетной программы и длительное в нем работать (задавать нагрузки, связи и т.д.), а встроить расчетную программу в нормальный графический редактор через автоматизацию невозможно).

В связи с вышеизложенным, назрел вопрос:
Как технологически наиболее правильно можно организовать разработку расчетного программного обеспечения с открытым исходным кодом?

Для совместной разработки кода создано общее хранилище на GitHub, используя которое каждый может поучаствовать в разработке :
https://github.com/chaosEagleOwl/source

На данным момент работа находится в стадии тестирования возможности совместной разработки.
Требования к программному обеспечению изложены в файле (ссылка README.md на GitHub): https://github.com/chaosEagleOwl/source/README.md

ТЗ на модуль формирования КЭ-сеток сформировано и помещено на GitHub.
На весь комплекс ТЗ формировать долго, видимо, будет чуть позже.

Сформирована доска для управления проектом, туда добавлены наиболее актуальные задачи.
Задачи проекта.

Последний раз редактировалось nickname2019, 06.10.2021 в 09:07.
Просмотров: 32831
 
Непрочитано 13.10.2021, 07:30
#301
Нубий-IV

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


Под VS2008/nanoCAD собирается, но требует более старую версию eigen - 3.1 (в текущей 3.4 не опознает новомодные шаблоны с переменным числом параметров, что, видимо, не лечится).
Похоже, VS2008 не опознает UTF-8 без BOM.
Миниатюры
Нажмите на изображение для увеличения
Название: 01.png
Просмотров: 314
Размер:	96.6 Кб
ID:	241689  
__________________
Учись долго, умри дураком

Последний раз редактировалось Нубий-IV, 13.10.2021 в 07:54.
Нубий-IV вне форума  
 
Непрочитано 13.10.2021, 08:19
#302
румата


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


Цитата:
Сообщение от nickname2019 Посмотреть сообщение
Может быть, работать с eigen через "обертки", на случай, если придется менять математическую библиотеку?
Абсолютно не понятно каким образом могут помочь обертки eigen при смене его на другую библиотеку.
румата вне форума  
 
Автор темы   Непрочитано 13.10.2021, 09:21
#303
nickname2019


 
Регистрация: 18.11.2019
Сообщений: 1,037


Цитата:
Сообщение от румата Посмотреть сообщение
Абсолютно не понятно каким образом могут помочь обертки eigen при смене его на другую библиотеку.
Можно доступ к функциям из eigen делать через прокладки:

код вызывает функцию-прокладку-> прокладка вызывает функцию eigen.

Тогда при смене математического модуля придется переписать только "прокладки", а основной код не трогать.
Иначе светит глобальная переработка, например у меня в процедурном типе tLocalStiffnessMatrixProcedure используется входной параметр MatrixXd & localStiffnessMatrix, что привязывает дальнейший код к eigen

typedef bool (*tLocalStiffnessMatrixProcedure)(const tFiniteElement & finiteElement, const tPoint3dArray& point3dArray, const tMaterialArray & materialArray, const tSectionArray &sectionArray, MatrixXd & localStiffnessMatrix);

Возможно, имеет смысл вводить промежуточный класс, который будет внутри себя использовать классы eigen.

----- добавлено через ~2 ч. -----
Цитата:
Сообщение от Нубий-IV Посмотреть сообщение
Под VS2008/nanoCAD собирается, но требует более старую версию eigen - 3.1 (в текущей 3.4 не опознает новомодные шаблоны с переменным числом параметров, что, видимо, не лечится).
Похоже, VS2008 не опознает UTF-8 без BOM.
Вы можете найти простой выход из ситуации или будем назад менять кодировку в основном файле? Хотелось бы, конечно, максимальной совместимости в проекте, чтобы везде собирался без проблем и кракозябр.

Последний раз редактировалось nickname2019, 13.10.2021 в 10:55.
nickname2019 вне форума  
 
Непрочитано 13.10.2021, 11:27
#304
румата


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


Цитата:
Сообщение от nickname2019 Посмотреть сообщение
Тогда при смене математического модуля придется переписать только "прокладки", а основной код не трогать.
Думается замена коней на переправе не очень хорошая затея. Что "прокладки" переписывать, что основной код - хрен редьки не слаще. СтОит все ж окончательно определится на основе какой математической библиотеки(или сразу двух таких либ) будем строить свой решатель. Как по мне, то для .Net-реализации решения и хранения матриц стОит использовать Math.NET (https://numerics.mathdotnet.com/) т.к. эта библиотека уже является оберткой над самым производительным ядром для решения задач линейной алгебры IntelMKL (https://numerics.mathdotnet.com/MKL.html). Проив eigen (BLAS) ничего не имею, но его будет нужно "оборачивать" самостоятельно, чего делать вообще не хочется.
румата вне форума  
 
Автор темы   Непрочитано 13.10.2021, 12:28
#305
nickname2019


 
Регистрация: 18.11.2019
Сообщений: 1,037


Цитата:
Сообщение от румата Посмотреть сообщение
Думается замена коней на переправе не очень хорошая затея. Что "прокладки" переписывать, что основной код - хрен редьки не слаще. СтОит все ж окончательно определится на основе какой математической библиотеки(или сразу двух таких либ) будем строить свой решатель. Как по мне, то для .Net-реализации решения и хранения матриц стОит использовать Math.NET (https://numerics.mathdotnet.com/) т.к. эта библиотека уже является оберткой над самым производительным ядром для решения задач линейной алгебры IntelMKL (https://numerics.mathdotnet.com/MKL.html). Проив eigen (BLAS) ничего не имею, но его будет нужно "оборачивать" самостоятельно, чего делать вообще не хочется.
Вроде бы MKL доступна бесплатно только для студентов и некоторых исследовательских задач (https://softline.ru/uploads/ad/0d/e5.../3a/origin.pdf).

Вот сравнение eigen с какой-то boost
https://journal.tusur.ru/storage/495...pdf?1472809983

Видимо, надо найти нормальную библиотеку с линейной алгеброй под C#.

Под с++, видимо, придется оставить eigen, но чтобы сделать код максимально совместимым с .Net (а также с целью не привязываться к конкретной библиотеке), классы Eigen я уберу из основного кода.

В коде будет некоторый (может быть абстрактный) класс типа tLinearAlgebra, который будет вызывать методы из самой математической библиотеки.

Последний раз редактировалось nickname2019, 13.10.2021 в 13:14.
nickname2019 вне форума  
 
Непрочитано 13.10.2021, 12:33
#306
Бахил

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


Offtop: О, уже до "решателя" добрались.
__________________
В конструктивных дискуссиях каждый участник укрепляется в своих заблуждениях.
Бахил на форуме  
 
Автор темы   Непрочитано 13.10.2021, 12:52
#307
nickname2019


 
Регистрация: 18.11.2019
Сообщений: 1,037


Цитата:
Сообщение от Бахил Посмотреть сообщение
Offtop: О, уже до "решателя" добрались.
По моему мнению на данный момент, объем требуемого кода - до 2 тыс. строк. Остальное - существующие библиотеки.
Может быть больше с учетом функций проверки и подбора сечений (там сейчас трудно объем оценить).
nickname2019 вне форума  
 
Непрочитано 13.10.2021, 14:52
#308
румата


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


Цитата:
Сообщение от nickname2019 Посмотреть сообщение
Вроде бы MKL доступна бесплатно только для студентов и некоторых исследовательских задач
Оно бесплатно для пользования(не коммерческого) в Windows, т.к. у создателей Math.NET есть на него лицензия и потому MKL идет в составе Math.NET.
Цитата:
Be aware that unlike the core of Math.NET Numerics including the native wrapper, which are both open source under the terms of the MIT/X11 license, the Intel MKL binaries themselves are closed source and non-free.

The Math.NET Numerics project does own an Intel MKL license (for Windows, no longer for Linux) and thus does have the right to distribute it along Math.NET Numerics. You can therefore use the Math.NET Numerics MKL native provider for free for your own use. However, it does not give you any right to redistribute it again yourself to customers of your own product. If you need to redistribute, buy a license from Intel. If unsure, contact the Intel sales team to clarify.
Перевод

Цитата:
Имейте в виду, что в отличие от ядра Math.NET Numerics, включая собственную оболочку, которые являются открытыми исходными кодами в соответствии с условиями лицензии MIT/X11, сами двоичные файлы Intel MKL являются закрытыми и несвободными.

Проект Math.NET Numerics имеет лицензию Intel MKL (для Windows, а не для Linux) и, таким образом, имеет право распространять ее вместе Math.NET Numerics. Таким образом, вы можете использовать Math.NET Numerics MKL native provider бесплатно для собственного использования. Однако это не дает вам никакого права перераспределять его снова самостоятельно среди клиентов собственного продукта. Если вам необходимо провести повторное распространение, приобретите лицензию у корпорации Intel. Если вы не уверены, обратитесь в отдел продаж Intel для уточнения.
румата вне форума  
 
Автор темы   Непрочитано 13.10.2021, 15:14
#309
nickname2019


 
Регистрация: 18.11.2019
Сообщений: 1,037


Цитата:
Сообщение от румата Посмотреть сообщение
Оно бесплатно для пользования(не коммерческого) в Windows, т.к. у создателей Math.NET есть на него лицензия и потому MKL идет в составе Math.NET.
Имхо, мы не можем его включать в свободный продукт, так как распространять нельзя. Наш код может предполагать любое, в т.ч. коммерческое использование.
nickname2019 вне форума  
 
Непрочитано 13.10.2021, 15:20
#310
румата


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


Цитата:
Сообщение от nickname2019 Посмотреть сообщение
Имхо, мы не можем его включать в свободный продукт
А мы его и не будем включать туда. Кому нужно сам его скачает и подключит для ускорения счета "адских" схем.
Цитата:
Сообщение от nickname2019 Посмотреть сообщение
Кроме того, наш код может предполагать любое, в т.ч. коммерческое использование.
Просто для коммерческого использования вместе с MKL нужно будет купить лицензию на MKL и все. Не так уж она дорого стоит.
румата вне форума  
 
Непрочитано 13.10.2021, 15:43
#311
Нубий-IV

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


Цитата:
Сообщение от nickname2019 Посмотреть сообщение
простой выход из ситуации
Мне помогло пересохранить как UTF с маркером BOM:
Если верить гуглу, разные версии компиляторов этот UTF - и BOMнутый, и неBOMнутый - ухитряются понимать по-разному. Вот VS2008 показывает файл без маркера правильно, а компилирует - неправильно; с маркером - все работает как надо. Другие студии - надо тестировать. (Пример - https://www.codeatcpp.com/2016/03/utf-8.html)
Цитата:
Сообщение от nickname2019 Посмотреть сообщение
максимальной совместимости в проекте, чтобы везде собирался без проблем
Как минимум условную компиляцию на разные кады или версии придется прописывать - выше уже приводил ошибки при сборке. Но это, наверное, по ходу дела будем добавлять, сейчас никто полный список не предскажет.

Будут как минимум проблемы при компиляции с eigen. Если для старых студий нужна старая версия, то в файле Types.h указывать полный путь к заголовкам (#include <C:\ArxLib\eigen\Eigen\Dense>) не годится: чтобы пересобраться под другую студию, надо будет каждый раз удалять старую папку с одной версией библиотеки и записывать другую версию.

Вообще предлагаю отвязать проекты от конкретных путей через свойства:
Фокус в том, что файл настроек окружения - это отдельный файл Environment.vsprops, в котором задана только пара макросов с путями, и каждому нужно поправить их под свое расположение библиотек; остальные настройки задаются по-прежнему в проекте, и их менять не придется. А пути к библиотекам после этого можно указывать просто - как и к другим стандартным "#include <Dense>".

Не соображу только, как настроить Git, чтобы он этот файл игнорировал: если включить его в gitignore сразу, то он не попадает в репозиторий, и его никто не скачает, чтобы настроить под себя; а если включить после добавления в репозиторий, то он уже не игнорируется. Как сделать, чтобы от с репозитория клонировался, но обратно после настройки уже не просился?
Миниатюры
Нажмите на изображение для увеличения
Название: 01.png
Просмотров: 281
Размер:	114.6 Кб
ID:	241701  Нажмите на изображение для увеличения
Название: 02.png
Просмотров: 281
Размер:	171.9 Кб
ID:	241704  
__________________
Учись долго, умри дураком
Нубий-IV вне форума  
 
Непрочитано 13.10.2021, 15:54
#312
румата


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


Цитата:
Сообщение от nickname2019 Посмотреть сообщение
Имхо, мы не можем его включать в свободный продукт, так как распространять нельзя.
Да, я тут еще вспомнил, что Math.NET, кроме собственного встроенного, использует еще и OpenBLAS в качестве нативного ядра линейной алгеьры, как альтернативу IntelMKL. Оно считает медленне MKL, но значительно быстрее встроенного.
румата вне форума  
 
Автор темы   Непрочитано 13.10.2021, 16:22
#313
nickname2019


 
Регистрация: 18.11.2019
Сообщений: 1,037


Проект, кажется, выходит на финальный этап, когда уже можно начинать писать по-делу.
Нужно обсудить, как в него вносить изменения. Видимо, делать merge через GitHubDesktop (как я все время делал) - не самая удачная идея, так как будет выбивать почву у людей, которые начали писать код на основе предыдущей версии.
После того, как румата поменял кодировку arxentrypoint, у меня синхронизация перестала закачиваться на gitHub (наверно, это нормально, так как кодировка файлов другая совершенно).
nickname2019 вне форума  
 
Непрочитано 13.10.2021, 16:49
#314
Нубий-IV

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


Руководство по Git: https://habr.com/ru/post/150673/ - коротенько, на 300 страниц. В пятой главе вроде про совместную работу пишут. Я туда через пару месяцев дочитаю, наверное - тогда тоже смогу поучаствовать.
__________________
Учись долго, умри дураком
Нубий-IV вне форума  
 
Непрочитано 13.10.2021, 16:59
| 1 #315
veb86

Проектировщик электрических сетей
 
Регистрация: 17.01.2014
Пенза
Сообщений: 197


Товарищи почитайте со стороны что вы пишите!
Румата не обижайся, тебе сказали С++, что ты к ним пристал со своим С#. Делай рядом свой проект, используй какие хочешь библиотеки. Я думаю ребята помогут тебе подключится к своему. Если у тебя будет что то интересное, то подсмотрят. У вас есть лидер, смерись с eigen.
Нубий-IV, вот у руматы C#, а у тебя бесплатный нанокад...
Я кончено не специалист по гиту, но с ZCAD-ом все просто, перед коммитом делаешь ПУЛЛ, решаешь все конфликты на этом моменте, делаешь коммит(помечаешь те файлы которые отправляешь) и пушешь (отправляешь). Естественно работать надо через тортоисеГИТ
Подробнее только подскажет zamtmn
veb86 вне форума  
 
Непрочитано 13.10.2021, 17:40
#316
румата


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


Цитата:
Сообщение от veb86 Посмотреть сообщение
Румата не обижайся, тебе сказали С++, что ты к ним пристал со своим С#.
Наше дело, как говорится, предложить. Ни к кому я не приставал. Я просто предложил готовое рабочее решение СЛАУ. Без танцев с бубнами вокруг сборки eigen.
Цитата:
Сообщение от veb86 Посмотреть сообщение
У вас есть лидер, смерись с eigen.
Не вопрос, пусть будет eigen. Главное что б кто-то дал ему лад.
Цитата:
Сообщение от nickname2019 Посмотреть сообщение
После того, как румата поменял кодировку arxentrypoint, у меня синхронизация перестала закачиваться на gitHub (наверно, это нормально, так как кодировка файлов другая совершенно).
Работать с репозиторием лучше всего прямо из VS. Лучше всего создать свою ветку кода и "ляпать" в ней. Все коммиты должны синхронизироваться независимо от кодировки файлов. А если коммитов не делать, то ничего и не будет синхронизироваться.
румата вне форума  
 
Непрочитано 13.10.2021, 18:03
| 1 #317
zamtmn

КИПиА
 
Регистрация: 21.03.2005
Tyumen
Сообщений: 1,497
<phrase 1=


>>Проект, кажется, выходит на финальный этап
оптимистично)))
zamtmn вне форума  
 
Непрочитано 13.10.2021, 19:41
#318
Сергей812


 
Регистрация: 10.08.2013
Сообщений: 9,463


Offtop:
Цитата:
Сообщение от nickname2019 Посмотреть сообщение
Проект, кажется, выходит на финальный этап, когда уже можно начинать писать по-делу.
ТЗ начинает приобретать какие то видимые черты скорее уж...
Сергей812 вне форума  
 
Непрочитано 15.10.2021, 09:02
#319
Бахил

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


Какое ТЗ? Полноценный проект. Осталось найти подрядчика. Можно и "хозспособом" реализовать.
__________________
В конструктивных дискуссиях каждый участник укрепляется в своих заблуждениях.
Бахил на форуме  
 
Непрочитано 24.10.2021, 08:25
2 | 1 #320
Нубий-IV

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


Чуток поигрался с этой ЭйГейской библиотекой.
Код:
[Выделить все]
 
// Тест скорости вычисления матрицы жесткости треугольной балки-стенки
#include <iostream>
#include <ctime>
#include <Dense>

using namespace std;

// Обозначения по книге
// МКЭ в проектировании транспортных сооружений.1981.Городецкий
// Стр.34
double E = 3e7;
double mu = 0.2;
double a = 1.1;
double b = 1.2;
double c = 1.3;

Eigen::MatrixXd K(6,6);

// Тестирующая функция
typedef void (*test_function)();
void test(test_function f)
{
    unsigned start_time = clock();
    for(int i = 0; i < 1000000; i++)
        f();
    unsigned end_time = clock();
    cout << "Time: " << (end_time - start_time) << " msec\n";
}


// В лоб по формулам из книги без оптимизации
void f1()
{
    double k = E / (1 - mu*mu);

    K(0,0) = k*(((2*a-b)*b-a*a)*mu+2*c*c+b*(b-2*a)+a*a)/(4*a*c);
    K(0,1) = k*((a-b)*mu-b+a)/(4*a);
    K(0,2) = k*(b*(b-a)*mu-2*c*c+(a-b)*b)/(4*a*c);
    K(0,3) = k*((b+a)*mu+b-a)/(4*a);
    K(0,4) = k*((a-b)*mu+b-a)/(4*c);
    K(0,5) = -k*mu/2;

    K(1,0) = K(1,0);
    K(1,1) = k*(-c*c*mu+c*c+b*(2*b-4*a)+2*a*a)/(4*a*c);
    K(1,2) = k*((b-2*a)*mu+b)/(4*a);
    K(1,3) = k*(c*c*mu-c*c+(2*a-2*b)*b)/(4*a*c);
    K(1,4) = k*(mu-1)/4;
    K(1,5) = k*(b-a)/(2*c);

    K(2,0) = K(0,2);
    K(2,1) = K(1,2);
    K(2,2) = k*(-b*b*mu+2*c*c+b*b)/(4*a*c);
    K(2,3) = k*(-b*mu-b)/(4*a);
    K(2,4) = k*(b*mu-b)/(4*c);
    K(2,5) = k*mu/2;

    K(3,0) = K(0,3);
    K(3,1) = K(1,3);
    K(3,2) = K(2,3);
    K(3,3) = k*(-c*c*mu+c*c+2*b*b)/(4*a*c);
    K(3,4) = k*(1-mu)/4;
    K(3,5) = -k*b/(2*c);

    K(4,0) = K(0,4);
    K(4,1) = K(1,4);
    K(4,2) = K(2,4);
    K(4,3) = K(3,4);
    K(4,4) = k*(a-a*mu)/(4*c);
    K(4,5) = 0;

    K(5,0) = K(0,5);
    K(5,1) = K(1,5);
    K(5,2) = K(2,5);
    K(5,3) = K(3,5);
    K(5,4) = K(4,5);
    K(5,5) = k*a/(2*c);
}

// В формулы, полученные и оптимизированные в Maxima
void f2()
{
    double k = E / (1 - mu*mu);
    double t1 = 1/a;
    double t2 = 1/c;
    double t3 = a*a;
    double t4 = 2*a*b;
    double t5 = b*b;
    double t6 = -t5;
    double t7 = c*c;
    double t8 = -2*t7;
    double t9 = t5*mu;
    double t10 = -a;
    double t11 = b+t10;
    double t12 = mu+1;
    double t13 = -(t1*t11*t12)/4;
    double t14 = (t1*t2*(t9-a*b*mu+t8+t6+a*b))/4;
    double t15 = b*mu;
    double t16 = (t1*(t15+a*mu+b+t10))/4;
    double t17 = mu-1;
    double t18 = -(t11*t2*t17)/4;
    double t19 = -mu/2;
    double t20 = -2*t5;
    double t21 = -t7;
    double t22 = t7*mu;
    double t23 = (t1*(t15-2*a*mu+b))/4;
    double t24 = (t1*t2*(t22+t21+t20+t4))/4;
    double t25 = t17/4;
    double t26 = (t11*t2)/2;
    double t27 = -(t1*b*t12)/4;
    double t28 = (b*t2*t17)/4;
    double t29 = mu/2;
    double t30 = -t17/4;
    double t31 = -(b*t2)/2;

    K(0,0) = k*(-(t1*t2*(t9-2*a*b*mu+t3*mu+t8+t6+t4-t3))/4);
    K(0,1) = k*(t13);
    K(0,2) = k*(t14);
    K(0,3) = k*(t16);
    K(0,4) = k*(t18);
    K(0,5) = k*(t19);

    K(1,0) = K(1,0);
    K(1,1) = k*(-(t1*t2*(t22+t21+t20+4*a*b-2*t3))/4);
    K(1,2) = k*(t23);
    K(1,3) = k*(t24);
    K(1,4) = k*(t25);
    K(1,5) = k*(t26);

    K(2,0) = K(0,2);
    K(2,1) = K(1,2);
    K(2,2) = k*(-(t1*t2*(t9+t8+t6))/4);
    K(2,3) = k*(t27);
    K(2,4) = k*(t28);
    K(2,5) = k*(t29);

    K(3,0) = K(0,3);
    K(3,1) = K(1,3);
    K(3,2) = K(2,3);
    K(3,3) = k*(-(t1*t2*(t22+t21+t20))/4);
    K(3,4) = k*(t30);
    K(3,5) = k*(t31);

    K(4,0) = K(0,4);
    K(4,1) = K(1,4);
    K(4,2) = K(2,4);
    K(4,3) = K(3,4);
    K(4,4) = k*(-(a*t2*t17)/4);
    K(4,5) = 0;

    K(5,0) = K(0,5);
    K(5,1) = K(1,5);
    K(5,2) = K(2,5);
    K(5,3) = K(3,5);
    K(5,4) = K(4,5);
    K(5,5) = k*((a*t2)/2);
}

// Численное интегрирование по одной точке
void f3()
{
    double u = 1.0/3;
    double w = 1.0/3;
    double W = 0.5;

    double x = a*u + b*w;
    double y = c*w;

    Eigen::MatrixXd InvFuxy(6,6);
    InvFuxy(0,0) = 1;
    InvFuxy(0,1) = 0;
    InvFuxy(0,2) = 0;
    InvFuxy(0,3) = 0;
    InvFuxy(0,4) = 0;
    InvFuxy(0,5) = 0;

    InvFuxy(1,0) = -1/a;
    InvFuxy(1,1) = 0;
    InvFuxy(1,2) = 1/a;
    InvFuxy(1,3) = 0;
    InvFuxy(1,4) = 0;
    InvFuxy(1,5) = 0;

    InvFuxy(2,0) = (b-a)/(a*c);
    InvFuxy(2,1) = 0;
    InvFuxy(2,2) = -b/(a*c);
    InvFuxy(2,3) = 0;
    InvFuxy(2,4) = 1/c;
    InvFuxy(2,5) = 0;

    InvFuxy(3,0) = 0;
    InvFuxy(3,1) = 1;
    InvFuxy(3,2) = 0;
    InvFuxy(3,3) = 0;
    InvFuxy(3,4) = 0;
    InvFuxy(3,5) = 0;

    InvFuxy(4,0) = 0;
    InvFuxy(4,1) = -1/a;
    InvFuxy(4,2) = 0;
    InvFuxy(4,3) = 1/a;
    InvFuxy(4,4) = 0;
    InvFuxy(4,5) = 0;

    InvFuxy(5,0) = 0;
    InvFuxy(5,1) = (b-a)/(a*c);
    InvFuxy(5,2) = 0;
    InvFuxy(5,3) = -b/(a*c);
    InvFuxy(5,4) = 0;
    InvFuxy(5,5) = 1/c;

    Eigen::MatrixXd dfxy(3,6);
    dfxy(0,0) = 0;
    dfxy(0,1) = 1;
    dfxy(0,2) = 0;
    dfxy(0,3) = 0;
    dfxy(0,4) = 0;
    dfxy(0,5) = 0;

    dfxy(1,0) = 0;
    dfxy(1,1) = 0;
    dfxy(1,2) = 0;
    dfxy(1,3) = 0;
    dfxy(1,4) = 0;
    dfxy(1,5) = 1;

    dfxy(2,0) = 0;
    dfxy(2,1) = 0;
    dfxy(2,2) = 1;
    dfxy(2,3) = 0;
    dfxy(2,4) = 1;
    dfxy(2,5) = 0;

    Eigen::MatrixXd fexy = dfxy * InvFuxy;

    double k = E / (1 - mu*mu);
    Eigen::MatrixXd D(3,3);
    D(0,0) = k;
    D(0,1) = k*mu;
    D(0,2) = 0;

    D(1,0) = k*mu;
    D(1,1) = k;
    D(1,2) = 0;

    D(2,0) = 0;
    D(2,1) = 0;
    D(2,2) = k*(1-mu)/2;

    double DetJ = a*c;

    K = fexy.transpose() * D * fexy * DetJ * W;
}

// Численное интегрирование - оптимизировано с учетом нулей. Сразу вычисляется fexy
void f4()
{
    double u = 1.0/3;
    double w = 1.0/3;
    double W = 0.5;

    double x = a*u + b*w;
    double y = c*w;

    Eigen::MatrixXd fexy(3,6);
    fexy(0,0) = -1/a;
    fexy(0,1) = 0;
    fexy(0,2) = 1/a;
    fexy(0,3) = 0;
    fexy(0,4) = 0;
    fexy(0,5) = 0;

    fexy(1,0) = 0;
    fexy(1,1) = (b-a)/(a*c);
    fexy(1,2) = 0;
    fexy(1,3) = -b/(a*c);
    fexy(1,4) = 0;
    fexy(1,5) = 1/c;

    fexy(2,0) = (b-a)/(a*c);
    fexy(2,1) = -1/a;
    fexy(2,2) = -b/(a*c);
    fexy(2,3) = 1/a;
    fexy(2,4) = 1/c;
    fexy(2,5) = 0;

    double k = E / (1 - mu*mu);
    Eigen::MatrixXd D(3,3);
    D(0,0) = k;
    D(0,1) = k*mu;
    D(0,2) = 0;

    D(1,0) = k*mu;
    D(1,1) = k;
    D(1,2) = 0;

    D(2,0) = 0;
    D(2,1) = 0;
    D(2,2) = k*(1-mu)/2;

    double DetJ = a*c;

    K = fexy.transpose() * D * fexy * DetJ * W;
}

// Численное интегрирование. Массивы C-style. Перемножение вручную с учетом нулей.
double R[6][6];
void cout_R()
{
    cout << "K =\n";
    for(int i = 0; i < 6; i++) {
        for(int j = 0; j < 6; j++)
            cout << '\t' << R[i][j];
        cout << '\n';
    }
}

void f5()
{
    double u = 1.0/3;
    double w = 1.0/3;
    double W = 0.5;

    double x = a*u + b*w;
    double y = c*w;

    double ac = a*c;
    double ba = b-a;

    double fexy02 = 1/a;
    double fexy00 = -fexy02;
    double fexy11 = ba/ac;
    double fexy13 = -fexy11;
    double fexy15 = 1/c;
    double fexy20 = fexy11; // оптимизировать: не хранить повторы
    double fexy21 = fexy00;
    double fexy22 = fexy13;
    double fexy23 = fexy02;
    double fexy24 = fexy15;

    double k = E / (1 - mu*mu);
    double D00 = k;
    double D01 = k*mu;
    double D10 = D01;
    double D11 = D00;
    double D22 = k*(1-mu)/2;

    double Dfexy00 = D00*fexy00;
    double Dfexy01 = D01*fexy11;
    double Dfexy02 = D00*fexy02;
    double Dfexy03 = D01*fexy13;
    double Dfexy05 = D01*fexy15;

    double Dfexy11 = D11*fexy11;
    double Dfexy12 = D10*fexy02;
    double Dfexy13 = D11*fexy13;
    double Dfexy15 = D11*fexy15;

    double Dfexy20 = D22*fexy20;
    double Dfexy21 = D22*fexy21;
    double Dfexy22 = D22*fexy22;
    double Dfexy23 = D22*fexy23;
    double Dfexy24 = D22*fexy24;

    double DetJ = a*c;
    double WDetJ = W*DetJ;

    R[0][0] = WDetJ*(Dfexy20*fexy20+Dfexy00*fexy00);
    R[0][1] = WDetJ*(Dfexy21*fexy20+Dfexy01*fexy00);
    R[0][2] = WDetJ*(Dfexy22*fexy20+Dfexy02*fexy00);
    R[0][3] = WDetJ*(Dfexy23*fexy20+Dfexy03*fexy00);
    R[0][4] = WDetJ*(Dfexy24*fexy20);
    R[0][5] = WDetJ*(Dfexy05*fexy00);

    R[1][0] = R[0][1];
    R[1][1] = WDetJ*(Dfexy21*fexy21+Dfexy11*fexy11);
    R[1][2] = WDetJ*(Dfexy22*fexy21+Dfexy12*fexy11);
    R[1][3] = WDetJ*(Dfexy23*fexy21+Dfexy13*fexy11);
    R[1][4] = WDetJ*(Dfexy24*fexy21);
    R[1][5] = WDetJ*(Dfexy15*fexy11);

    R[2][0] = R[0][2];
    R[2][1] = R[1][2];
    R[2][2] = WDetJ*(Dfexy22*fexy22+Dfexy02*fexy02);
    R[2][3] = WDetJ*(Dfexy23*fexy22+Dfexy03*fexy02);
    R[2][4] = WDetJ*(Dfexy24*fexy22);
    R[2][5] = WDetJ*(Dfexy05*fexy02);

    R[3][0] = R[0][3];
    R[3][1] = R[1][3];
    R[3][2] = R[2][3];
    R[3][3] = WDetJ*(Dfexy23*fexy23+Dfexy13*fexy13);
    R[3][4] = WDetJ*(Dfexy24*fexy23);
    R[3][5] = WDetJ*(Dfexy15*fexy13);

    R[4][0] = R[0][4];
    R[4][1] = R[1][4];
    R[4][2] = R[2][4];
    R[4][3] = R[3][4];
    R[4][4] = WDetJ*(Dfexy24*fexy24);
    R[4][5] = 0;

    R[5][0] = R[0][5];
    R[5][1] = R[1][5];
    R[5][2] = R[2][5];
    R[5][3] = R[3][5];
    R[5][4] = R[4][5];
    R[5][5] = WDetJ*(Dfexy15*fexy15);
}

int main()
{
    cout << "Function test\n";
    f1();
    cout << "K = \n" << K << endl;

    f2();
    cout << "K = \n" << K << endl;

    f3();
    cout << "K = \n" << K << endl;

    f4();
    cout << "K = \n" << K << endl;

    f5();
    cout_R();

    cout << "Speed test\n";
    test(f1);
    test(f2);
    test(f3);
    test(f4);
    test(f5);

    return 0;
}
Результаты. Числитель - VS2008 + Eigen-3.1 / Знаменатель - GCC8 + Eigen-3.4.
  • 16/47
    Аналитическое решение.
    Готовые формулы для элементов матрицы, взятые из книги. Наверное, можно еще чуть ускорить, если вынести повторы и сгруппировать умножения.
  • 16/47
    Аналитическое решение.
    Те же формулы, полученные в символьном виде в Maxima, после функции Optimize(). Можно еще ускорить, если присваивать одноразовые временные переменные сразу в матрицу.
  • 1015/797
    Численное интегрирование.
    Eigen - реализация формул МКЭ "Как в книге" - вычисление векторов координат, потом векторов коэффициентов, потом матрицы градиентов, и наконец матрицы жесткости [ K ]=[ B^T ][ D ][ B ].
  • 625/391
    Численное интегрирование.
    Eigen - попытка оптимизации, за счет избавления от перемножения матриц с большим числом нулей. Так же сделано в примере на хабре из поста 137 - матрица градиентов вычисляется сразу, а не перемножением других матриц.
  • 15/16
    Численное интегрирование.
    C-массивы. Ручное перемножение матриц с учетом нулей, симметрии и без лишнего прожорливого транспонирования.
Вопросы:
  • Результаты зависят от компилятора. Может быть, другие Студии покажут другие результаты. Как бы не пришлось еще под каждый компилятор версии писать.
  • Забавно, что ручное численное интегрирование обставило аналитику. Может, в более сложных элементах (меньше нулей, больше точек интегрирования, переменные жескости) будет больший объем счета, но там и формулы будут сильно сложнее.
  • Типовая формула из учебников [ K ]=[ B^T ][ D ][ B ] что-то сильно тормозит. Даже с избыточным перемножением нулевых элементов, неучетом симметрии и лишним транспонированием - непонятно, откуда разница в десятки раз. Может, там Евгену что-то подкрутить надо в настройках?

Второй тест - решение систем.
Код:
[Выделить все]
 
// EigenTest.cpp : Defines the entry point for the console application.
//

#include "stdafx.h"


int _tmain(int argc, _TCHAR* argv[])
{
	int MatrixSize = 200000;
	int BandWidth = 50;

	std::cout << "Stiffness matrix calulation...\n";
	unsigned start_time = std::clock();

	Eigen::SparseMatrix<double> K(MatrixSize, MatrixSize);
	std::vector<Eigen::Triplet<double>> T;
	for(int i = 0; i < MatrixSize; i++) {
		int bandMinIndex = std::max(0, i - BandWidth + 1);
		int bandMaxIndex = std::min(i + BandWidth, MatrixSize);
		for(int j = bandMinIndex; j < bandMaxIndex; j++) {
			double distance = i-j;
			double max_distance = BandWidth;
			double relativeDistance = distance / max_distance;
			T.push_back(Eigen::Triplet<double>(i, j, 1.0 - std::pow(relativeDistance, 2)));
		}
	}
	K.setFromTriplets(T.begin(), T.end());

	unsigned end_time = std::clock();
	std::cout << (end_time - start_time) / 1000.0 << " seconds\n";
	std::cout << "...done.\n\n";
	//K.coeffRef(0,0) = 0.0;

	std::cout << "Forces vector calulation...\n";
	start_time = std::clock();
	Eigen::VectorXd F(MatrixSize);

	for(int i = 0; i < MatrixSize; i++)
		F[i] = (double)i / MatrixSize;

	end_time = std::clock();
	std::cout << (end_time - start_time) / 1000.0 << " seconds\n";
	std::cout << "...done.\n\n";

	std::cout << "Equation system solving...\n";
	start_time = std::clock();

	Eigen::SimplicialLDLT<Eigen::SparseMatrix<double>> solver(K);
	Eigen::VectorXd x = solver.solve(F);

	end_time = std::clock();
	std::cout << (end_time - start_time) / 1000.0 << " seconds\n";
	std::cout << "...done.\n";

	//std::cout << "Solver info: " << solver.info() << std::endl;
	//std::cout << "\nStiffness:\n" << K << std::endl;
	//std::cout << "\nForces:\n" << F << std::endl;
	//std::cout << "\nDisplacements:\n" << x << std::endl;

	return 0;
}
Есть вопросы:
  • Матрицы держит в памяти (вроде даже без учета симметрии - хранит и верх, и низ). Значит, большие схемы на слабом железе считать не выйдет?
  • Не понимает, что ноль на диагонали - это автоматическая связь, а не ошибка расчета. Можно ему объяснить, что он неправ?
  • Приостановка/возобновление расчета невозможны?
__________________
Учись долго, умри дураком
Нубий-IV вне форума  
Ответ
Вернуться   Форум DWG.RU > Программное обеспечение > Программирование > Как создать расчетное программное обеспечение с открытым исходным кодом (конструктивные решения)

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

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


Похожие темы
Тема Автор Раздел Ответов Последнее сообщение
СП 335.1325800.2017 «Крупнопанельные конструктивные системы. Правила проектирования» (Обсуждение) Armin Прочее. Архитектура и строительство 37 07.11.2018 06:55
Фирменные решения по пропуску коммуникаций через стены подвала Regby Конструкции зданий и сооружений 2 07.04.2010 20:43
устройство и возможные конструктивные решения вентфасада из кирпича Ivansobaka Каменные и армокаменные конструкции 1 16.12.2009 06:38
Конструктивные решения по перемычкам в многослойных кирпичных стенах! Westroy Архитектура 16 30.11.2009 13:57
Конструктивные решения монтажных соединений многоэтажных зданий на высокопрочных болтах VoRoNoFF Конструкции зданий и сооружений 1 04.04.2009 00:41