[nickname2019]

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

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

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

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

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

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

[Нубий-IV]

Под VS2008/nanoCAD собирается, но требует более старую версию eigen - 3.1 (в текущей 3.4 не опознает новомодные шаблоны с переменным числом параметров, что, видимо, не лечится).

Похоже, VS2008 не опознает UTF-8 без BOM.

[румата]

Сообщение от nickname2019 Может быть, работать с eigen через "обертки", на случай, если придется менять математическую библиотеку?

Абсолютно не понятно каким образом могут помочь обертки eigen при смене его на другую библиотеку.

[nickname2019]

Сообщение от румата Абсолютно не понятно каким образом могут помочь обертки 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 Тогда при смене математического модуля придется переписать только "прокладки", а основной код не трогать.

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

[nickname2019]

Сообщение от румата Думается замена коней на переправе не очень хорошая затея. Что "прокладки" переписывать, что основной код - хрен редьки не слаще. СтОит все ж окончательно определится на основе какой математической библиотеки(или сразу двух таких либ) будем строить свой решатель. Как по мне, то для .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, который будет вызывать методы из самой математической библиотеки.

[Бахил]

Offtop: О, уже до "решателя" добрались.

[nickname2019]

Сообщение от Бахил Offtop: О, уже до "решателя" добрались.

По моему мнению на данный момент, объем требуемого кода - до 2 тыс. строк. Остальное - существующие библиотеки.
Может быть больше с учетом функций проверки и подбора сечений (там сейчас трудно объем оценить).

[румата]

Сообщение от 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 для уточнения.

[nickname2019]

Сообщение от румата Оно бесплатно для пользования(не коммерческого) в Windows, т.к. у создателей Math.NET есть на него лицензия и потому MKL идет в составе Math.NET.

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

[румата]

Сообщение от nickname2019 Имхо, мы не можем его включать в свободный продукт

А мы его и не будем включать туда. Кому нужно сам его скачает и подключит для ускорения счета "адских" схем.

Сообщение от nickname2019 Кроме того, наш код может предполагать любое, в т.ч. коммерческое использование.

Просто для коммерческого использования вместе с MKL нужно будет купить лицензию на MKL и все. Не так уж она дорого стоит.

[Нубий-IV]

Сообщение от 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 сразу, то он не попадает в репозиторий, и его никто не скачает, чтобы настроить под себя; а если включить после добавления в репозиторий, то он уже не игнорируется. Как сделать, чтобы от с репозитория клонировался, но обратно после настройки уже не просился?

[румата]

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

Да, я тут еще вспомнил, что Math.NET, кроме собственного встроенного, использует еще и OpenBLAS в качестве нативного ядра линейной алгеьры, как альтернативу IntelMKL. Оно считает медленне MKL, но значительно быстрее встроенного.

[nickname2019]

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

[Нубий-IV]

Руководство по Git: https://habr.com/ru/post/150673/ - коротенько, на 300 страниц. В пятой главе вроде про совместную работу пишут. Я туда через пару месяцев дочитаю, наверное - тогда тоже смогу поучаствовать.

[veb86]

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

[румата]

Сообщение от veb86 Румата не обижайся, тебе сказали С++, что ты к ним пристал со своим С#.

Наше дело, как говорится, предложить. Ни к кому я не приставал. Я просто предложил готовое рабочее решение СЛАУ. Без танцев с бубнами вокруг сборки eigen.

Сообщение от veb86 У вас есть лидер, смерись с eigen.

Не вопрос, пусть будет eigen. Главное что б кто-то дал ему лад.

Сообщение от nickname2019 После того, как румата поменял кодировку arxentrypoint, у меня синхронизация перестала закачиваться на gitHub (наверно, это нормально, так как кодировка файлов другая совершенно).

Работать с репозиторием лучше всего прямо из VS. Лучше всего создать свою ветку кода и "ляпать" в ней. Все коммиты должны синхронизироваться независимо от кодировки файлов. А если коммитов не делать, то ничего и не будет синхронизироваться.

[zamtmn]

>>Проект, кажется, выходит на финальный этап
оптимистично)))

[Сергей812]

Offtop:

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

ТЗ начинает приобретать какие то видимые черты скорее уж...

[Бахил]

Какое ТЗ? Полноценный проект. Осталось найти подрядчика. Можно и "хозспособом" реализовать.

[Нубий-IV]

Чуток поигрался с этой ЭйГейской библиотекой.

Код:

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

 
// Тест скорости вычисления матрицы жесткости треугольной балки-стенки
#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;
}

Есть вопросы:

• Матрицы держит в памяти (вроде даже без учета симметрии - хранит и верх, и низ). Значит, большие схемы на слабом железе считать не выйдет?
• Не понимает, что ноль на диагонали - это автоматическая связь, а не ошибка расчета. Можно ему объяснить, что он неправ?
• Приостановка/возобновление расчета невозможны?