[nickname2019]

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

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

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

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

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

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

[румата]

Сообщение от nickname2019 Если юзер забыл приложить связи на систему, то в процессе решения может быть ошибка (если матрицу решать гауссом, то на диагонали в процессе решения появиться ноль).

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

[nickname2019]

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

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

[Нубий-IV]

У меня тут неприятность образовалась в одном месте. Я неожиданно вспомнил, что нанокад - 32-битный. Получается, я сообственную свечку в собственный нанокад не смогу вставить - вот где геморрой-то будет.

Еще вопрос. Если мы собираемся узлы с разным числом переменных делать, то нельзя сразу сказать, какому индексу в матрице системы соответствует данная степень свободы данного КЭ. Надо сначала пройтись по всем КЭ и составить список использованных направлений, а для этого надо все их матрицы жесткости посчитать и в глобальные координаты перевести - только там все нули видны будут. И только после опроса последнего элемента можно записать в общую матрицу первый. И тогда либо все матрицы хранить надо, что приводит как минимум к утроению расхода памяти (и даже больше), потому что они а) не разреженные б) с нижним треугольником в) не записаны друг поверх друга в общих узлах г) висят в памяти параллельно со сборкой общей матрицы. Либо при сборке общей матрицы все элементы придется обсчитывать повторно. Либо надо у узла список направлений указывать заранее (а несовпадение при расчете считать ошибкой) - тогда в схеме надо еще и узлы как-то задавать, а не только элементы. Либо таки Лироскады правы, и в узлах просто 6 неизвестных, которые при расчете автозакрепляются, если надо.

Кстати, кто в курсе - в кадовскую базу элементов можно в многопоточном режиме обращаться, или это еще одно узкое место получается?

[Бахил]

Мгновенно-изменяемая определяется "на берегу" нулями на диагонали.
Вот вырожденную матрицу (без опорных закреплений) выявить труднее.
И, кстати, матрица без опорных закреплений не всегда вырожденная.

Сообщение от Нубий-IV все матрицы хранить надо

Они нужны и для определения внутренних усилий/напряжений. Можно хранить только ненулевые элементы.

[nickname2019]

Сообщение от Нубий-IV У меня тут неприятность образовалась в одном месте. Я неожиданно вспомнил, что нанокад - 32-битный. Получается, я сообственную свечку в собственный нанокад не смогу вставить - вот где геморрой-то будет.

Сольвер можно сделать отдельным exe, так как подключить 64-разрядную DLL к 32-разрядной системе вряд ли выйдет. Что делать с 64-разрядным акадом пока не понятно. Не понятно нужно ли создавать многопоточную задачу внутри его основного процесса.

Сообщение от Нубий-IV Еще вопрос. Если мы собираемся узлы с разным числом переменных делать, то нельзя сразу сказать, какому индексу в матрице системы соответствует данная степень свободы данного КЭ. Надо сначала пройтись по всем КЭ и составить список использованных направлений, а для этого надо все их матрицы жесткости посчитать и в глобальные координаты перевести - только там все нули видны будут. И только после опроса последнего элемента можно записать в общую матрицу первый. И тогда либо все матрицы хранить надо, что приводит как минимум к утроению расхода памяти (и даже больше), потому что они а) не разреженные б) с нижним треугольником в) не записаны друг поверх друга в общих узлах г) висят в памяти параллельно со сборкой общей матрицы.

Я тоже об этом думал. На следующей неделе попробую набросать код для расчета стержневых систем, матрицу жесткости для стержней я вроде бы нашел (см. вложение).
Сначала нужно определить тип элемента в зависимости от имен степеней свободы, потом формировать типы узлов в зависимости от примыкаемых элементов (до формирования матрицы).

Я думал матрицы жесткости формировать отдельно по каждым уровням (этажам), обращать, а потом как суперэлементы объединять. На многоэтажках с типовыми этажами должен быть серьезный эффект.

Сообщение от Нубий-IV Кстати, кто в курсе - в кадовскую базу элементов можно в многопоточном режиме обращаться, или это еще одно узкое место получается?

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

[Бахил]

Сообщение от nickname2019 (см. вложение).

МКЭ там и близко не стояло. Обычный метод перемещений.
МКЭ для стержня хорошо описан в Клаф, Пензиен. Динамика сооружений.
Ну а матрицы для пластин смотрите в Зенкевиче.
https://dwg.ru/dnl/4923 и https://dwg.ru/dnl/1949

[Нубий-IV]

Еще один тест скорости, более реалистичный, с элементом балки-стенки. Программа читает входной файл in.txt и записывает выходной out.pos. Формат входного файла - аналогичный статье из поста 137. Чтобы не мучиться с ручным вводом - дальше простейший скриптовый генератор. Код генератора сохранить с расширением HTA и запускать как обычную программу. Результаты смотреть в GMSH.

Код:

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

 
#include<iostream>
#include<fstream>
#include<vector>
#include<ctime>
#include<Dense>
#include<Sparse>


namespace Timer {
	unsigned start_time;
	void Start(const char* msg)  { std::cout << msg; start_time = std::clock(); }
	void Stop()  { std::cout << (std::clock() - start_time) / 1000.0 << " seconds\n"; }
}


struct Material {
	double mu;
	double E;
};


struct Section {
	double h;
};


struct Node {
	double x;
	double y;
	double z;
};


struct Element{
	int materialId;
	int sectionId;
	int nodesIds[3];
};


struct Constraint {
	int nodeId;
	int direction;
};


struct Load {
	int nodeId;
	int direction;
	double value;
};


typedef std::vector<Material> Materials;
typedef std::vector<Section> Sections;
typedef std::vector<Node> Nodes;
typedef std::vector<Element> Elements;
typedef std::vector<Constraint> Constraints;
typedef std::vector<Load> Loads;

typedef Eigen::Matrix<double, 6, 9> ElementMatrixT;
typedef Eigen::Matrix<double, 6, 6> ElementLocalMatrixK;
typedef Eigen::Matrix<double, 9, 9> ElementGlobalMatrixK;

typedef std::vector<ElementMatrixT *> ElementMatricesT;
typedef std::vector<ElementLocalMatrixK *> ElementLocalMatricesK;
typedef std::vector<ElementGlobalMatrixK *> ElementGlobalMatricesK;

typedef Eigen::SparseMatrix<double> SystemMatrixK;
typedef Eigen::VectorXd LoadsVector;
typedef Eigen::VectorXd Variables;

typedef Eigen::SimplicialLDLT<Eigen::SparseMatrix<double>> Solver;


void ReadMaterials(std::ifstream& infile, Materials& materials)
{
	int nMaterials;
	infile >> nMaterials;
	materials.resize(nMaterials);
	for(int i =0; i < nMaterials; i++)
		infile >> materials[i].mu >> materials[i].E;
}


void ReadSections(std::ifstream& infile, Sections& sections)
{
	int nSections;
	infile >> nSections;
	sections.resize(nSections);
	for(int i =0; i < nSections; i++)
		infile >> sections[i].h;
}


void ReadNodes(std::ifstream& infile, Nodes& nodes)
{
	int nNodes;
	infile >> nNodes;
	nodes.resize(nNodes);
	for(int i = 0; i < nNodes; i++)
		infile >> nodes[i].x >> nodes[i].y >> nodes[i].z;

}


void ReadElements(std::ifstream& infile, Elements& elements)
{
	int nElements;
	infile >> nElements;
	elements.resize(nElements);
	for(int i = 0; i < nElements; i++)
		infile >> elements[i].materialId >> elements[i].sectionId >> elements[i].nodesIds[0] >> elements[i].nodesIds[1] >> elements[i].nodesIds[2];
}


void ReadConstraints(std::ifstream& infile, Constraints& constraints)
{
	int nConstraints;
	infile >> nConstraints;
	constraints.resize(nConstraints);
	for(int i = 0; i < nConstraints; i++)
		infile >> constraints[i].nodeId >> constraints[i].direction;
}


void ReadLoads(std::ifstream& infile, Loads& loads)
{
	int nLoads;
	infile >> nLoads;
	loads.resize(nLoads);
	for(int i = 0; i < nLoads; i++)
		infile >> loads[i].nodeId >> loads[i].direction >> loads[i].value;
}


void CalcElementMatrixT(const Element& element, const Nodes& nodes, ElementMatrixT& t)
{
	Eigen::Vector3d x;
	x[0] = nodes[element.nodesIds[1]].x - nodes[element.nodesIds[0]].x;
	x[1] = nodes[element.nodesIds[1]].y - nodes[element.nodesIds[0]].y;
	x[2] = nodes[element.nodesIds[1]].z - nodes[element.nodesIds[0]].z;

	Eigen::Vector3d v;
	v[0] = nodes[element.nodesIds[2]].x - nodes[element.nodesIds[0]].x;
	v[1] = nodes[element.nodesIds[2]].y - nodes[element.nodesIds[0]].y;
	v[2] = nodes[element.nodesIds[2]].z - nodes[element.nodesIds[0]].z;

	Eigen::Vector3d z = x.cross(v);
	Eigen::Vector3d y = z.cross(x);

	x.normalize();
	y.normalize();

	t(0,0) = x[0]; t(0,1) = x[1]; t(0,2) = x[2]; t(0,3) = 0;    t(0,4) =  0;   t(0,5) = 0;    t(0,6) = 0;    t(0,7) = 0;    t(0,8) = 0;
	t(1,0) = y[0]; t(1,1) = y[1]; t(1,2) = y[2]; t(1,3) = 0;    t(1,4) =  0;   t(1,5) = 0;    t(1,6) = 0;    t(1,7) = 0;    t(1,8) = 0;
	t(2,0) = 0;    t(2,1) = 0;    t(2,2) = 0;    t(2,3) = x[0]; t(2,4) = x[1]; t(2,5) = x[2]; t(2,6) = 0;    t(2,7) = 0;    t(2,8) = 0;
	t(3,0) = 0;    t(3,1) = 0;    t(3,2) = 0;    t(3,3) = y[0]; t(3,4) = y[1]; t(3,5) = y[2]; t(3,6) = 0;    t(3,7) = 0;    t(3,8) = 0;
	t(4,0) = 0;    t(4,1) = 0;    t(4,2) = 0;    t(4,3) = 0;    t(4,4) = 0;    t(4,5) = 0;    t(4,6) = x[0]; t(4,7) = x[1]; t(4,8) = x[2];
	t(5,0) = 0;    t(5,1) = 0;    t(5,2) = 0;    t(5,3) = 0;    t(5,4) = 0;    t(5,5) = 0;    t(5,6) = y[0]; t(5,7) = y[1]; t(5,8) = y[2];
}

void CalcElementMatricesT(const Elements& elements, const Nodes& nodes, ElementMatricesT& matricesT)
{
	matricesT.resize(elements.size());

	for(unsigned i = 0; i < elements.size(); i++) {
		ElementMatrixT* t = new ElementMatrixT;
		CalcElementMatrixT(elements[i], nodes, *t);
		matricesT[i] = t;
	}
}


void CalcElementLocalMatrixK(const Element& element, const Materials& materials, const Sections& sections, const Nodes& nodes, const ElementMatrixT& t, ElementLocalMatrixK& kl)
{
	Eigen::VectorXd globalCoord(9);

	globalCoord[0] = 0;
	globalCoord[1] = 0;
	globalCoord[2] = 0;

	globalCoord[3] = nodes[element.nodesIds[1]].x - nodes[element.nodesIds[0]].x;
	globalCoord[4] = nodes[element.nodesIds[1]].y - nodes[element.nodesIds[0]].y;
	globalCoord[5] = nodes[element.nodesIds[1]].z - nodes[element.nodesIds[0]].z;

	globalCoord[6] = nodes[element.nodesIds[2]].x - nodes[element.nodesIds[0]].x;
	globalCoord[7] = nodes[element.nodesIds[2]].y - nodes[element.nodesIds[0]].y;
	globalCoord[8] = nodes[element.nodesIds[2]].z - nodes[element.nodesIds[0]].z;

	Eigen::VectorXd localCoord = t * globalCoord;

	double a = localCoord[2]; // x2
	double b = localCoord[4]; // x3
	double c = localCoord[5]; // y3

	double E = materials[element.materialId].E;
	double mu = materials[element.materialId].mu;
	double h = sections[element.sectionId].h;
	double f = E * h / (1 - mu*mu);

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

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

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

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

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

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

void CalcElementLocalMatricesK(const Elements& elements, const Materials& materials, const Sections& sections, const Nodes& nodes, 
			       const ElementMatricesT& elementMatricesT, ElementLocalMatricesK& elementLocalMatricesK)
{
	elementLocalMatricesK.resize(elements.size());

	for(unsigned i = 0; i < elements.size(); i++) {
		ElementLocalMatrixK* kl = new ElementLocalMatrixK;
		CalcElementLocalMatrixK(elements[i], materials, sections, nodes, *elementMatricesT[i], *kl);
		elementLocalMatricesK[i] = kl;
	}
	
}


void CalcElementGlobalMatrixK(const ElementLocalMatrixK& kl, const ElementMatrixT& t, ElementGlobalMatrixK& kg)
{
	kg = t.transpose() * kl * t;
}

void CalcElementGlobalMatricesK(const ElementMatricesT& elementMatricesT, const ElementLocalMatricesK& elementLocalMatricesK, ElementGlobalMatricesK& elementGlobalMatricesK)
{
	elementGlobalMatricesK.resize(elementLocalMatricesK.size());

	for(unsigned i = 0; i < elementLocalMatricesK.size(); i++) {
		ElementGlobalMatrixK* elementGlobalMatrixK = new ElementGlobalMatrixK;
		CalcElementGlobalMatrixK(*elementLocalMatricesK[i], *elementMatricesT[i], *elementGlobalMatrixK);
		elementGlobalMatricesK[i] = elementGlobalMatrixK;
	}
}


void PlaceElementMatrixK(const Element& element, const ElementGlobalMatrixK& k, SystemMatrixK& K)
{
	for(int i = 0; i < k.rows(); i++) {
		int iNodeId = i / 3;
		int iDirectionId = i % 3;
		for(int j = 0; j < k.cols(); j++) {
			int jNodeId = j / 3;
			int jDirectionId = j % 3;
			int m = 6 * element.nodesIds[iNodeId] + iDirectionId;
			int n = 6 * element.nodesIds[jNodeId] + jDirectionId;
			K.coeffRef(m,n) += k(i,j);
		}
	}
}


void CalcSystemMatrixK(const Elements& elements, const ElementGlobalMatricesK&  elementGlobalMatricesK, SystemMatrixK& K)
{
	for(unsigned i = 0; i < elements.size(); i++) {
		PlaceElementMatrixK(elements[i], *elementGlobalMatricesK[i], K);
	}
}


void SetConstraint(const Constraint constraint, SystemMatrixK& K)
{
	K.coeffRef(6 * constraint.nodeId + constraint.direction,6 * constraint.nodeId + constraint.direction) = 1e10;
}


void CalcConstraints(const Constraints& constraints, SystemMatrixK& K)
{
	for(unsigned i = 0; i < constraints.size(); i++)
		SetConstraint(constraints[i], K);
}


void SetLoad(const Load load, LoadsVector& L)
{
	L.coeffRef(6 * load.nodeId + load.direction) = load.value;
}


void CalcLoads(const Loads& loads, LoadsVector& L)
{
	L.setZero();
	for(unsigned i = 0; i < loads.size(); i++)
		SetLoad(loads[i], L);
}


void WriteU(std::ofstream& outfile, const Nodes& nodes, const Elements& elements, const Variables& U)
{
	outfile << "View \"Ux\"{\n";
	for(unsigned i = 0; i < elements.size(); i++)
		outfile << "ST(" <<
			nodes[elements[i].nodesIds[0]].x << ',' << nodes[elements[i].nodesIds[0]].y << ',' << nodes[elements[i].nodesIds[0]].z << ',' <<
			nodes[elements[i].nodesIds[1]].x << ',' << nodes[elements[i].nodesIds[1]].y << ',' << nodes[elements[i].nodesIds[1]].z << ',' <<
			nodes[elements[i].nodesIds[2]].x << ',' << nodes[elements[i].nodesIds[2]].y << ',' << nodes[elements[i].nodesIds[2]].z << "){" <<
			U[6*elements[i].nodesIds[0]] << "," << U[6*elements[i].nodesIds[1]] << "," << U[6*elements[i].nodesIds[2]] << "};\n";
	outfile << "};\n";

	outfile << "View \"Uy\"{\n";
	for(unsigned i = 0; i < elements.size(); i++)
		outfile << "ST(" <<
			nodes[elements[i].nodesIds[0]].x << ',' << nodes[elements[i].nodesIds[0]].y << ',' << nodes[elements[i].nodesIds[0]].z << ',' <<
			nodes[elements[i].nodesIds[1]].x << ',' << nodes[elements[i].nodesIds[1]].y << ',' << nodes[elements[i].nodesIds[1]].z << ',' <<
			nodes[elements[i].nodesIds[2]].x << ',' << nodes[elements[i].nodesIds[2]].y << ',' << nodes[elements[i].nodesIds[2]].z << "){" <<
			U[6*elements[i].nodesIds[0]+1] << "," << U[6*elements[i].nodesIds[1]+1] << "," << U[6*elements[i].nodesIds[2]+1] << "};\n";
	outfile << "};\n";

	outfile << "View \"Uz\"{\n";
	for(unsigned i = 0; i < elements.size(); i++)
		outfile << "ST(" <<
			nodes[elements[i].nodesIds[0]].x << ',' << nodes[elements[i].nodesIds[0]].y << ',' << nodes[elements[i].nodesIds[0]].z << ',' <<
			nodes[elements[i].nodesIds[1]].x << ',' << nodes[elements[i].nodesIds[1]].y << ',' << nodes[elements[i].nodesIds[1]].z << ',' <<
			nodes[elements[i].nodesIds[2]].x << ',' << nodes[elements[i].nodesIds[2]].y << ',' << nodes[elements[i].nodesIds[2]].z << "){" <<
			U[6*elements[i].nodesIds[0]+2] << "," << U[6*elements[i].nodesIds[1]+2] << "," << U[6*elements[i].nodesIds[2]+2] << "};\n";
	outfile << "};\n";
}


int main(int argc, char* argv[])
{
	std::ifstream infile("in.txt");
	std::ofstream outfile("out.pos");

	Timer::Start("Read data...");
		Materials materials;
		ReadMaterials(infile, materials);

		Sections sections;
		ReadSections(infile, sections);

		Nodes nodes;
		ReadNodes(infile, nodes);

		Elements elements;
		ReadElements(infile, elements);

		Constraints constraints;
		ReadConstraints(infile, constraints);

		Loads loads;
		ReadLoads(infile, loads);
	Timer::Stop();

	Timer::Start("Element transformation matrices calculation...");
		ElementMatricesT elementMatricesT;
		CalcElementMatricesT(elements, nodes, elementMatricesT);
	Timer::Stop();

	Timer::Start("Element local stiffness matrices calculation...");
		ElementLocalMatricesK elementLocalMatricesK;
		CalcElementLocalMatricesK(elements, materials, sections, nodes, elementMatricesT, elementLocalMatricesK);
	Timer::Stop();

	Timer::Start("Element global stiffness matrices calculation...");
		ElementGlobalMatricesK elementGlobalMatricesK;
		CalcElementGlobalMatricesK(elementMatricesT, elementLocalMatricesK, elementGlobalMatricesK);
	Timer::Stop();

	Timer::Start("System siffnes matrix calulation...");
		int nVars = 6 * nodes.size();
		SystemMatrixK K(nVars, nVars);
		CalcSystemMatrixK(elements, elementGlobalMatricesK, K);
	Timer::Stop();

	Timer::Start("Constraints calulation...");
		CalcConstraints(constraints, K);
	Timer::Stop();

	Timer::Start("Loads calulation...");
		LoadsVector L(nVars);
		CalcLoads(loads, L);
	Timer::Stop();

	Timer::Start("Equations system solving...");
		Solver solver(K);
		Variables U = solver.solve(L);
	Timer::Stop();

	Timer::Start("Results writting...");
		WriteU(outfile, nodes, elements, U);
	Timer::Stop();
}

Код:

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

<h1>Генератор исходных данных</h1>
<p>Ширина балки-стенки <input name="i_b"  value="1"> м</p>
<p>Высота балки-стенки <input name="i_h"  value="1"> м</p>
<p>Число элементов по ширине <input name="i_nb" value="1"></p>
<p>Число элементов по высоте <input name="i_nh" value="1"></p>
<p>Модуль упругости <input name="i_E"  value="3e7"> Кпа</p>
<p>Коэффициент Пуассона <input name="i_mu" value="0.2"></p>
<p>Толщина балки-стенки <input name="i_t"  value="0.5"> м</p>
<p>Распределенная нагрузка по верху балки-стенки <input name="i_q"  value="100"> Кн/м</p>
<input type="button" value="Генерировать" onclick="Generate()">

<script>

function Generate()
{
	fso = new ActiveXObject("Scripting.FileSystemObject")
	var file = fso.OpenTextFile("In.txt", 2, true)

	b  = parseFloat(i_b.value)
	nb = parseFloat(i_nb.value)
	h  = parseFloat(i_h.value)
	nh = parseFloat(i_nh.value)
	E  = parseFloat(i_E.value)
	mu = parseFloat(i_mu.value)
	t  = parseFloat(i_t.value)
	q  = parseFloat(i_q.value)

	// Материалы
	file.Write("1\n" + mu + " " + E + "\n")
	
	// Сечения
	file.Write("1\n" + t + "\n")
	
	// Узлы
	nNodes = (nb + 1) * (nh + 1)
	file.Write("" + nNodes + "\n")
	sizeX = b/nb
	sizeY = h/nh
	for(j = 0; j <= nh; j++)
		for(i = 0; i <= nb; i++) {
			x = i * sizeX
			y = j * sizeY
			z = 0
			file.Write("" + x + " " + y + " " + z + "\n")
		}

	// Элементы
	nElements = 2 * nb * nh
	file.Write("" + nElements + "\n")
	for(j = 0; j < nh; j++)
		for(i = 0; i < nb; i++) {
			node0 = i + (nb + 1) * j
			node1 = i + (nb + 1) * j + 1
			node2 = i + (nb + 1) * (j + 1)
			file.Write("0 0 " + node0 + " " + node1 + " " + node2 + "\n")

			node0 = i + (nb + 1) * (j + 1)
			node1 = i + (nb + 1) * j + 1
			node2 = i + (nb + 1) * (j + 1) + 1
			file.Write("0 0 " + node0 + " " + node1 + " " + node2 + "\n")
		}

	// Связи
	nConstraints = 6 + 5 * nb + 4 * (nb + 1) * nh
	file.Write("" + nConstraints + "\n")
	file.Write("0 0\n0 1\n0 2\n0 3\n0 4\n0 5\n")
	for(i = 1; i <= nb; i++) {
		file.Write("" + i + " 1\n")
		file.Write("" + i + " 2\n")
		file.Write("" + i + " 3\n")
		file.Write("" + i + " 4\n")
		file.Write("" + i + " 5\n")
	}

	for(i = nb + 1; i < (nb + 1) * (nh + 1); i++) {
		file.Write("" + i + " 2\n")
		file.Write("" + i + " 3\n")
		file.Write("" + i + " 4\n")
		file.Write("" + i + " 5\n")
	}

	// Нагрузки
	nLoads = nb + 1
	file.Write("" + nLoads + "\n")

	f1 = q * sizeX
	f2 = q * sizeX / 2
	
	node = (nb + 1) * nh
	file.Write("" + node + " 1 " + f2 + "\n")

	for(i = (nb + 1) * nh + 1; i < (nb + 1) * (nh + 1) - 1; i++)
		file.Write("" + i + " 1 " + f1 + "\n")	

	node = (nb + 1) * (nh + 1) - 1
	file.Write("" + node + " 1 " + f2 + "\n")

	file.Close()	
}
</script>

Однако, самая прожорливая операция - вставка в матрицу системы. Даже оптимизация "нули не вставлять" разгоняет ее вдвое. Это какой-то позор.

Сообщение от nickname2019 Сольвер можно сделать отдельным exe

Кстати, так он больше похож на ОбПеньСорс. Если у него будет интерфейс для ввода данных и чтения результатов:

• В виде функций c++
• В виде функций c#
• В виде функций Python
• В виде текстового языка
• В виде двоичного формата

То поучаствовать смогут те, кто заходил в тему, и советовал выше это все.

[румата]

Сообщение от Нубий-IV А у гибридного элемента?

Нашел у себя вот такое


----- добавлено через ~11 мин. -----

Сообщение от Нубий-IV Формат входного файла - аналогичный статье из поста 137.

Думается стОит принять XML-формат для файла входных данных по причине возможности его валидации XML -схемой.

[Нубий-IV]

Однако, гугл помогл. Оказывается, у разреженной матрицы есть метод reserve(), в который можно передать профиль матрицы - вектор, у которого для каждого столбца указано, сколько ненулевых элементов в столбце будет. Тогда память выделится только один раз. Я в запас просуммировал количество всех элементов на этапе сборки глобальных матриц элементов, не учитывая наложений. Где без резервирования было 12с - с резервированием стало 0.016с; ускорение в 750 раз.

Код:

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

 
#include<iostream>
#include<fstream>
#include<vector>
#include<ctime>
#include<Dense>
#include<Sparse>


namespace Timer {
	unsigned start_time;
	void Start(const char* msg)  { std::cout << msg; start_time = std::clock(); }
	void Stop()  { std::cout << (std::clock() - start_time) / 1000.0 << " seconds\n"; }
}


struct Material {
	double mu;
	double E;
};


struct Section {
	double h;
};


struct Node {
	double x;
	double y;
	double z;
};


struct Element{
	int materialId;
	int sectionId;
	int nodesIds[3];
};


struct Constraint {
	int nodeId;
	int direction;
};


struct Load {
	int nodeId;
	int direction;
	double value;
};


typedef std::vector<Material> Materials;
typedef std::vector<Section> Sections;
typedef std::vector<Node> Nodes;
typedef std::vector<Element> Elements;
typedef std::vector<Constraint> Constraints;
typedef std::vector<Load> Loads;

typedef Eigen::Matrix<double, 6, 9> ElementMatrixT;
typedef Eigen::Matrix<double, 6, 6> ElementLocalMatrixK;
typedef Eigen::Matrix<double, 9, 9> ElementGlobalMatrixK;

typedef std::vector<ElementMatrixT *> ElementMatricesT;
typedef std::vector<ElementLocalMatrixK *> ElementLocalMatricesK;
typedef std::vector<ElementGlobalMatrixK *> ElementGlobalMatricesK;

typedef Eigen::VectorXi SystemMatrixProfile;
typedef Eigen::SparseMatrix<double> SystemMatrixK;
typedef Eigen::VectorXd LoadsVector;
typedef Eigen::VectorXd Variables;

typedef Eigen::SimplicialLDLT<Eigen::SparseMatrix<double>> Solver;


void ReadMaterials(std::ifstream& infile, Materials& materials)
{
	int nMaterials;
	infile >> nMaterials;
	materials.resize(nMaterials);
	for(int i =0; i < nMaterials; i++)
		infile >> materials[i].mu >> materials[i].E;
}


void ReadSections(std::ifstream& infile, Sections& sections)
{
	int nSections;
	infile >> nSections;
	sections.resize(nSections);
	for(int i =0; i < nSections; i++)
		infile >> sections[i].h;
}


void ReadNodes(std::ifstream& infile, Nodes& nodes)
{
	int nNodes;
	infile >> nNodes;
	nodes.resize(nNodes);
	for(int i = 0; i < nNodes; i++)
		infile >> nodes[i].x >> nodes[i].y >> nodes[i].z;

}


void ReadElements(std::ifstream& infile, Elements& elements)
{
	int nElements;
	infile >> nElements;
	elements.resize(nElements);
	for(int i = 0; i < nElements; i++)
		infile >> elements[i].materialId >> elements[i].sectionId >> elements[i].nodesIds[0] >> elements[i].nodesIds[1] >> elements[i].nodesIds[2];
}


void ReadConstraints(std::ifstream& infile, Constraints& constraints)
{
	int nConstraints;
	infile >> nConstraints;
	constraints.resize(nConstraints);
	for(int i = 0; i < nConstraints; i++)
		infile >> constraints[i].nodeId >> constraints[i].direction;
}


void ReadLoads(std::ifstream& infile, Loads& loads)
{
	int nLoads;
	infile >> nLoads;
	loads.resize(nLoads);
	for(int i = 0; i < nLoads; i++)
		infile >> loads[i].nodeId >> loads[i].direction >> loads[i].value;
}


void CalcElementMatrixT(const Element& element, const Nodes& nodes, ElementMatrixT& t)
{
	Eigen::Vector3d x;
	x[0] = nodes[element.nodesIds[1]].x - nodes[element.nodesIds[0]].x;
	x[1] = nodes[element.nodesIds[1]].y - nodes[element.nodesIds[0]].y;
	x[2] = nodes[element.nodesIds[1]].z - nodes[element.nodesIds[0]].z;

	Eigen::Vector3d v;
	v[0] = nodes[element.nodesIds[2]].x - nodes[element.nodesIds[0]].x;
	v[1] = nodes[element.nodesIds[2]].y - nodes[element.nodesIds[0]].y;
	v[2] = nodes[element.nodesIds[2]].z - nodes[element.nodesIds[0]].z;

	Eigen::Vector3d z = x.cross(v);
	Eigen::Vector3d y = z.cross(x);

	x.normalize();
	y.normalize();

	t(0,0) = x[0]; t(0,1) = x[1]; t(0,2) = x[2]; t(0,3) = 0;    t(0,4) =  0;   t(0,5) = 0;    t(0,6) = 0;    t(0,7) = 0;    t(0,8) = 0;
	t(1,0) = y[0]; t(1,1) = y[1]; t(1,2) = y[2]; t(1,3) = 0;    t(1,4) =  0;   t(1,5) = 0;    t(1,6) = 0;    t(1,7) = 0;    t(1,8) = 0;
	t(2,0) = 0;    t(2,1) = 0;    t(2,2) = 0;    t(2,3) = x[0]; t(2,4) = x[1]; t(2,5) = x[2]; t(2,6) = 0;    t(2,7) = 0;    t(2,8) = 0;
	t(3,0) = 0;    t(3,1) = 0;    t(3,2) = 0;    t(3,3) = y[0]; t(3,4) = y[1]; t(3,5) = y[2]; t(3,6) = 0;    t(3,7) = 0;    t(3,8) = 0;
	t(4,0) = 0;    t(4,1) = 0;    t(4,2) = 0;    t(4,3) = 0;    t(4,4) = 0;    t(4,5) = 0;    t(4,6) = x[0]; t(4,7) = x[1]; t(4,8) = x[2];
	t(5,0) = 0;    t(5,1) = 0;    t(5,2) = 0;    t(5,3) = 0;    t(5,4) = 0;    t(5,5) = 0;    t(5,6) = y[0]; t(5,7) = y[1]; t(5,8) = y[2];
}


void CalcElementMatricesT(const Elements& elements, const Nodes& nodes, ElementMatricesT& matricesT)
{
	matricesT.resize(elements.size());

	for(unsigned i = 0; i < elements.size(); i++) {
		ElementMatrixT* t = new ElementMatrixT;
		CalcElementMatrixT(elements[i], nodes, *t);
		matricesT[i] = t;
	}
}


void CalcElementLocalMatrixK(const Element& element, const Materials& materials, const Sections& sections, const Nodes& nodes, const ElementMatrixT& t, ElementLocalMatrixK& kl)
{
	Eigen::VectorXd globalCoord(9);

	globalCoord[0] = 0;
	globalCoord[1] = 0;
	globalCoord[2] = 0;

	globalCoord[3] = nodes[element.nodesIds[1]].x - nodes[element.nodesIds[0]].x;
	globalCoord[4] = nodes[element.nodesIds[1]].y - nodes[element.nodesIds[0]].y;
	globalCoord[5] = nodes[element.nodesIds[1]].z - nodes[element.nodesIds[0]].z;

	globalCoord[6] = nodes[element.nodesIds[2]].x - nodes[element.nodesIds[0]].x;
	globalCoord[7] = nodes[element.nodesIds[2]].y - nodes[element.nodesIds[0]].y;
	globalCoord[8] = nodes[element.nodesIds[2]].z - nodes[element.nodesIds[0]].z;

	Eigen::VectorXd localCoord = t * globalCoord;

	double a = localCoord[2]; // x2
	double b = localCoord[4]; // x3
	double c = localCoord[5]; // y3

	double E = materials[element.materialId].E;
	double mu = materials[element.materialId].mu;
	double h = sections[element.sectionId].h;
	double f = E * h / (1 - mu*mu);

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

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

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

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

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

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


void CalcElementLocalMatricesK(const Elements& elements, const Materials& materials, const Sections& sections, const Nodes& nodes, 
			       const ElementMatricesT& elementMatricesT, ElementLocalMatricesK& elementLocalMatricesK)
{
	elementLocalMatricesK.resize(elements.size());

	for(unsigned i = 0; i < elements.size(); i++) {
		ElementLocalMatrixK* kl = new ElementLocalMatrixK;
		CalcElementLocalMatrixK(elements[i], materials, sections, nodes, *elementMatricesT[i], *kl);
		elementLocalMatricesK[i] = kl;
	}
	
}


void CalcElementGlobalMatrixK(const ElementLocalMatrixK& kl, const ElementMatrixT& t, ElementGlobalMatrixK& kg)
{
	kg = t.transpose() * kl * t;
}


void CalcElementGlobalMatricesK(const Elements& elements, const ElementMatricesT& elementMatricesT, const ElementLocalMatricesK& elementLocalMatricesK, ElementGlobalMatricesK& elementGlobalMatricesK, SystemMatrixProfile& profile)
{
	elementGlobalMatricesK.resize(elementLocalMatricesK.size());

	for(unsigned i = 0; i < elementLocalMatricesK.size(); i++) {
		ElementGlobalMatrixK* elementGlobalMatrixK = new ElementGlobalMatrixK;
		CalcElementGlobalMatrixK(*elementLocalMatricesK[i], *elementMatricesT[i], *elementGlobalMatrixK);
		elementGlobalMatricesK[i] = elementGlobalMatrixK;

		for(int j = 0; j < elementGlobalMatrixK->cols(); j++) {
			int jNodeId = j / 3;
			int jDirectionId = j % 3;
			int n = 6 * elements[i].nodesIds[jNodeId] + jDirectionId;
			profile[n] += elementGlobalMatrixK->rows();
		}
	}
}


void PlaceElementMatrixK(const Element& element, const ElementGlobalMatrixK& k, SystemMatrixK& K)
{
	for(int i = 0; i < k.rows(); i++) {
		int iNodeId = i / 3;
		int iDirectionId = i % 3;
		for(int j = 0; j < k.cols(); j++) {
			int jNodeId = j / 3;
			int jDirectionId = j % 3;
			int m = 6 * element.nodesIds[iNodeId] + iDirectionId;
			int n = 6 * element.nodesIds[jNodeId] + jDirectionId;
			K.coeffRef(m,n) += k(i,j);
		}
	}
}


void CalcSystemMatrixK(const Elements& elements, const ElementGlobalMatricesK&  elementGlobalMatricesK, SystemMatrixK& K)
{
	for(unsigned i = 0; i < elements.size(); i++) {
		PlaceElementMatrixK(elements[i], *elementGlobalMatricesK[i], K);
	}
}


void SetConstraint(const Constraint constraint, SystemMatrixK& K)
{
	K.coeffRef(6 * constraint.nodeId + constraint.direction,6 * constraint.nodeId + constraint.direction) = 1e10;
}


void CalcConstraints(const Constraints& constraints, SystemMatrixK& K)
{
	for(unsigned i = 0; i < constraints.size(); i++)
		SetConstraint(constraints[i], K);
}


void SetLoad(const Load load, LoadsVector& L)
{
	L.coeffRef(6 * load.nodeId + load.direction) = load.value;
}


void CalcLoads(const Loads& loads, LoadsVector& L)
{
	L.setZero();
	for(unsigned i = 0; i < loads.size(); i++)
		SetLoad(loads[i], L);
}


void WriteU(std::ofstream& outfile, const Nodes& nodes, const Elements& elements, const Variables& U)
{
	outfile << "View \"Ux\"{\n";
	for(unsigned i = 0; i < elements.size(); i++)
		outfile << "ST(" <<
			nodes[elements[i].nodesIds[0]].x << ',' << nodes[elements[i].nodesIds[0]].y << ',' << nodes[elements[i].nodesIds[0]].z << ',' <<
			nodes[elements[i].nodesIds[1]].x << ',' << nodes[elements[i].nodesIds[1]].y << ',' << nodes[elements[i].nodesIds[1]].z << ',' <<
			nodes[elements[i].nodesIds[2]].x << ',' << nodes[elements[i].nodesIds[2]].y << ',' << nodes[elements[i].nodesIds[2]].z << "){" <<
			U[6*elements[i].nodesIds[0]] << "," << U[6*elements[i].nodesIds[1]] << "," << U[6*elements[i].nodesIds[2]] << "};\n";
	outfile << "};\n";

	outfile << "View \"Uy\"{\n";
	for(unsigned i = 0; i < elements.size(); i++)
		outfile << "ST(" <<
			nodes[elements[i].nodesIds[0]].x << ',' << nodes[elements[i].nodesIds[0]].y << ',' << nodes[elements[i].nodesIds[0]].z << ',' <<
			nodes[elements[i].nodesIds[1]].x << ',' << nodes[elements[i].nodesIds[1]].y << ',' << nodes[elements[i].nodesIds[1]].z << ',' <<
			nodes[elements[i].nodesIds[2]].x << ',' << nodes[elements[i].nodesIds[2]].y << ',' << nodes[elements[i].nodesIds[2]].z << "){" <<
			U[6*elements[i].nodesIds[0]+1] << "," << U[6*elements[i].nodesIds[1]+1] << "," << U[6*elements[i].nodesIds[2]+1] << "};\n";
	outfile << "};\n";

	outfile << "View \"Uz\"{\n";
	for(unsigned i = 0; i < elements.size(); i++)
		outfile << "ST(" <<
			nodes[elements[i].nodesIds[0]].x << ',' << nodes[elements[i].nodesIds[0]].y << ',' << nodes[elements[i].nodesIds[0]].z << ',' <<
			nodes[elements[i].nodesIds[1]].x << ',' << nodes[elements[i].nodesIds[1]].y << ',' << nodes[elements[i].nodesIds[1]].z << ',' <<
			nodes[elements[i].nodesIds[2]].x << ',' << nodes[elements[i].nodesIds[2]].y << ',' << nodes[elements[i].nodesIds[2]].z << "){" <<
			U[6*elements[i].nodesIds[0]+2] << "," << U[6*elements[i].nodesIds[1]+2] << "," << U[6*elements[i].nodesIds[2]+2] << "};\n";
	outfile << "};\n";
}


int main(int argc, char* argv[])
{
	std::ifstream infile("in.txt");
	std::ofstream outfile("out.pos");

	Timer::Start("Read data...");
		Materials materials;
		ReadMaterials(infile, materials);

		Sections sections;
		ReadSections(infile, sections);

		Nodes nodes;
		ReadNodes(infile, nodes);

		Elements elements;
		ReadElements(infile, elements);

		Constraints constraints;
		ReadConstraints(infile, constraints);

		Loads loads;
		ReadLoads(infile, loads);
	Timer::Stop();

	Timer::Start("Element transformation matrices calculation...");
		ElementMatricesT elementMatricesT;
		CalcElementMatricesT(elements, nodes, elementMatricesT);
	Timer::Stop();

	Timer::Start("Element local stiffness matrices calculation...");
		ElementLocalMatricesK elementLocalMatricesK;
		CalcElementLocalMatricesK(elements, materials, sections, nodes, elementMatricesT, elementLocalMatricesK);
	Timer::Stop();

	Timer::Start("Element global stiffness matrices calculation...");
		int nVars = 6 * nodes.size();
		SystemMatrixProfile profile(nVars);
		for(int i = 0; i < nVars; i++)
			profile[i]=1;

		ElementGlobalMatricesK elementGlobalMatricesK;
		CalcElementGlobalMatricesK(elements, elementMatricesT, elementLocalMatricesK, elementGlobalMatricesK, profile);
	Timer::Stop();

	Timer::Start("System siffnes matrix calulation...");
		SystemMatrixK K(nVars, nVars);
		K.reserve(profile);
		CalcSystemMatrixK(elements, elementGlobalMatricesK, K);
	Timer::Stop();

	Timer::Start("Constraints calulation...");
		CalcConstraints(constraints, K);
	Timer::Stop();

	Timer::Start("Loads calulation...");
		LoadsVector L(nVars);
		CalcLoads(loads, L);
	Timer::Stop();

	Timer::Start("Equations system solving...");
		Solver solver(K);
		Variables U = solver.solve(L);
	Timer::Stop();

	Timer::Start("Results writting...");
		WriteU(outfile, nodes, elements, U);
	Timer::Stop();
}

А вот тест на 120000 узлов (как у свечки):

И это без оптимизаий - хранится верх и низ матрицы, добавляются в т.ч. нулевые элементы (хотя у оболочек нулей будет меньше). Таким уже можно пользоваться.

Сообщение от румата XML-формат для файла входных данных

Скорее уж, основным должен быть либо двоичный файл, либо програмный интерфейс. А к нему разные утилиты импорта - экспорта: хоть XML, хоть Лировский TXT, хоть HTML - страничка.

[румата]

Сообщение от Нубий-IV Скорее уж, основным должен быть либо двоичный файл, либо програмный интерфейс. А к нему разные утилиты импорта - экспорта: хоть XML, хоть Лировский TXT, хоть HTML - страничка.

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

[Нубий-IV]

Сообщение от румата Ввод данных через бинарный файл точно не нужно делать.

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

Код:

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

 
#include<iostream>
#include<fstream>
#include<ctime>
#include<windows.h>


namespace Timer {
	unsigned start_time;
	void Start(const char* msg)  { std::cout << msg << std::endl; start_time = std::clock(); }
	void Stop()  { std::cout << (std::clock() - start_time) / 1000.0 << " seconds\n\n"; }
}


struct Node {
	Node():x(1), y(2), z(3) {}
	double x;
	double y;
	double z;
};


void WriteFormatted(int nNodes, Node* nodes)
{
	std::ofstream outfile("File.txt");
	
	for(int i = 0; i < nNodes; i++)
		outfile << "Node  =  { " << nodes[i].x << " , " << nodes[i].y << " , " << nodes[i].z << " }\n";
}


void WriteBinary(int nNodes, Node* nodes)
{
	HANDLE hFile = CreateFile("File.bin", GENERIC_WRITE, 0, 0, CREATE_ALWAYS, FILE_ATTRIBUTE_NORMAL, 0);

	WriteFile(hFile, nodes, nNodes * sizeof(Node), 0, 0);

	CloseHandle(hFile);
}


void main()
{
	const int nNodes = 10000000;
	Node* nodes = new Node[nNodes];

	Timer::Start("WriteFormatted");
		WriteFormatted(nNodes, nodes);
	Timer::Stop();

	Timer::Start("WriteBinary");
		WriteBinary(nNodes, nodes);
	Timer::Stop();
}

При форматированном выводе получится, что сначала ARX-модуль 20с записывает данные, потом решатель 20с их читает, делает расчет за 30с (см. тест выше), 20с пишет результаты, и ARX их обратно 20с читает - получится 80с возни с файлами и 30с расчета. И, чем сложнее формат, тем больше потери на чтение-запись. А в двоичной форме - 10с обработки файлов и 30с расчета. Еще обмен данными между программами через двоичный файл подкачки поддерживается системой (CreateFileMapping). И никто не мешает наклепать сторонних конвертеров под любые форматы. Но выбирать медленный формат как основной - по мне, неудачное решение.

У меня вообще, как у непрограммиста, проблема в том, что я синтаксис языка немного понимаю, а как его правильно применять - не очень. Даже простейший вопрос "как организовать перебор точек Гаусса" вызывает сложности.

Код:

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

 
#include<iostream>
#include<ctime>


namespace Timer {
	unsigned start_time;
	void Start(const char* msg)  { std::cout << msg << std::endl; start_time = std::clock(); }
	void Stop()  { std::cout << (std::clock() - start_time) / 1000.0 << " seconds\n\n"; }
}


// Число повторов в цикле
const int N = 1000000000;


// Подынтегральная функция
double f(double u)
{
	return 0.5 * u;
}


// Точка Гаусса - число
// Простейшее решение. Эталон скорости
// При смене числа точек требует правки в двух местах: заголовок цикла и строка суммирования
namespace Test1 {
	void run() {
		double integral = 0;
		for(int i = 0; i < N; i++) {
			integral += 2.0 * f(0.0);
		}
		std::cout << "integral = " << integral << std::endl;
	}
}


// Точка Гаусса - число в массиве
// При нескольких точках придется использовать массивы и циклы
// При смене числа точек требует правки в двух местах: заголовок цикла и строка суммирования
namespace Test2 {
	void run() {
		const int nPoints = 1;
		double u[nPoints] = {0.0};
		double W[nPoints] = {2.0};

		double integral = 0;
		for(int i = 0; i < N; i++) {
			for(int n = 0; n < nPoints; n++) {
				integral += W[n] * f(u[n]);
			}
		}
		std::cout << "integral = " << integral << std::endl;
	}
}


// Точка Гаусса - статическая структура
// Группировка координат и весов в структуру
// При смене числа точек требует правки в двух местах: заголовок цикла и строка суммирования
namespace Test3 {
	struct GaussPoint {
		double u;
		double W;
	} gaussPoint = { 0.0, 2.0 };

	const int nPoints = 1;

	void run() {
		double integral = 0;
		for(int i = 0; i < N; i++) {
			for(int n = 0; n < nPoints; n++) {
				integral += gaussPoint.W * f(gaussPoint.u);
			}
		}
		std::cout << "integral = " << integral << std::endl;
	}
}


// Точка Гаусса - структура, возвращаемая функцией
// При смене числа точек требует правки в двух местах: заголовок цикла и строка суммирования
namespace Test4 {
	struct GaussPoint {
		double u;
		double W;
	};

	GaussPoint GetPoint()
	{
		GaussPoint p = {0.0, 2.0};
		return p;
	}

	const int nPoints = 1;

	void run() {
		double integral = 0;
		for(int i = 0; i < N; i++) {
			for(int n = 0; n < nPoints; n++) {
				GaussPoint p = GetPoint();
				integral += p.W * f(p.u);
			}
		}
		std::cout << "integral = " << integral << std::endl;
	}
}


// Точка Гаусса - структура, возвращаемая функцией из массива
// При нескольких точках нужны массивы точек
// При смене числа точек требует правки в двух местах: заголовок цикла и строка суммирования
namespace Test5 {
	struct GaussPoint {
		double u;
		double W;
	};

	const int nPoints = 1;

	GaussPoint GetPoint()
	{
		static GaussPoint points[nPoints] = {0.0, 2.0};
		return points[0];
	}

	void run() {
		double integral = 0;
		for(int i = 0; i < N; i++) {
			for(int n = 0; n < nPoints; n++) {
				GaussPoint p = GetPoint();
				integral += p.W * f(p.u);
			}
		}
		std::cout << "integral = " << integral << std::endl;
	}
}


// Точка Гаусса - струкура с конструктором
// При смене числа точек требует правки в двух местах: заголовок цикла и строка суммирования
namespace Test6 {
	struct GaussPoint {
		GaussPoint(double u, double W): u(u), W(W) {}
		double u;
		double W;
	};

	const int nPoints = 1;

	void run() {
		double integral = 0;
		for(int i = 0; i < N; i++) {
			for(int n = 0; n < nPoints; n++) {
				GaussPoint p(0.0, 2.0);
				integral += p.W * f(p.u);
			}
		}
		std::cout << "integral = " << integral << std::endl;
	}
}


// Точка Гаусса - класс
// Вместо структур - классы. Вместо прямого доступа к полям - функции
// При смене числа точек требует правки в двух местах: заголовок цикла и строка суммирования
namespace Test7 {
	class GaussPoint {
	public:
		GaussPoint(double u, double W): m_u(u), m_W(W) {}
		double u() { return m_u; }
		double W() { return m_W; }
	private:
		double m_u;
		double m_W;
	};

	const int nPoints = 1;

	void run() {
		double integral = 0;
		for(int i = 0; i < N; i++) {
			for(int n = 0; n < nPoints; n++) {
				GaussPoint p(0.0, 2.0);
				integral += p.W() * f(p.u());
			}
		}
		std::cout << "integral = " << integral << std::endl;
	}
}


// Точка Гаусса - класс с приватными полями + класс схемы интегрирования
// Точки защищены от случайных изменений
// При смене числа точек требует единственного исправления - строки получения схемы
namespace Test8 {
	class GaussPoint {
	public:
		GaussPoint(double u, double W): m_u(u), m_W(W) {}
		double u() { return m_u; }
		double W() { return m_W; }
	private:
		double m_u;
		double m_W;
	};

	class GaussScheme {
	public:
		int NumPoints() { return 1; }
		GaussPoint Point() { return m_gaussPoint; }
	private:
		static GaussPoint m_gaussPoint;
	};

	GaussPoint GaussScheme::m_gaussPoint = GaussPoint(0.0, 2.0); 

	void run() {
		double integral = 0;
		for(int i = 0; i < N; i++) {
			GaussScheme scheme;
			for(int n = 0; n < scheme.NumPoints(); n++) {
				GaussPoint p = scheme.Point();
				integral += p.W() * f(p.u());
			}
		}
		std::cout << "integral = " << integral << std::endl;
	}
}


// Точка Гаусса - класс в массиве + класс схемы интегрирования
// При нескольких точках нужен массив
// При смене числа точек требует единственного исправления - строки получения схемы
namespace Test9 {
	class GaussPoint {
	public:
		GaussPoint(double u, double W): m_u(u), m_W(W) {}
		double u() { return m_u; }
		double W() { return m_W; }
	private:
		double m_u;
		double m_W;
	};

	class GaussScheme {
	public:
		int NumPoints() { return 1; }
		GaussPoint Point(int n) { return m_gaussPoints[n]; }
	private:
		static GaussPoint m_gaussPoints[];
	};

	GaussPoint GaussScheme::m_gaussPoints[] = {GaussPoint(0.0, 2.0)}; 

	void run() {
		double integral = 0;
		for(int i = 0; i < N; i++) {
			GaussScheme scheme;
			for(int n = 0; n < scheme.NumPoints(); n++) {
				GaussPoint p = scheme.Point(n);
				integral += p.W() * f(p.u());
			}
		}
		std::cout << "integral = " << integral << std::endl;
	}
}


// Точка Гаусса и схема интегрирования - абстрактные классы
// При смене числа точек требует единственного исправления - строки получения схемы
// Тормозное решение.
namespace Test10 {
	namespace Gauss {
		class Point {
		public:
			Point(double u, double W): m_u(u), m_W(W) {}
			double u() { return m_u; }
			double W() { return m_W; }
		private:
			double m_u;
			double m_W;
		};

		class IScheme {
		public:
			virtual int NumPoints() = 0;
			virtual Point* GetPoint(int n) = 0;
		};

		class Scheme: public IScheme {
		public:
			virtual int NumPoints() { return 1; }
			virtual Point* GetPoint(int n) { return &points[0]; }
		private:
			static Point points[1];
		} scheme;

		Point Scheme::points[1] = { Point(0.0, 2.0) };

		IScheme* GetScheme()
		{
			return &scheme;
		}
	}

	void run() {
		double integral = 0;
		for(int i = 0; i < N; i++) {
			Gauss::IScheme* scheme = Gauss::GetScheme();
			for(int n = 0; n < scheme->NumPoints(); n++) {
				Gauss::Point* point = scheme->GetPoint(n);
				integral += point->W() * f(point->u());
			}
		}
		std::cout << "integral = " << integral << std::endl;
	}
}


// Точка Гаусса - структура. Схема интегрирования - пространство имен.
// При смене числа точек требует правки в двух местах: заголовок цикла и строка суммирования
namespace Test11 {
	namespace Gauss {
		struct Point1D {
			double u;
			double W;
		};

		namespace Line1P {
			const int Points = 1;

			Point1D Point(int n)
			{
				static Point1D points[Points] = {0.0, 2.0};
				return points[n];
			}

		}
	}

	void run() {
		double integral = 0;
		for(int i = 0; i < N; i++) {
			for(int n = 0; n < Gauss::Line1P::Points; n++) {
				Gauss::Point1D point = Gauss::Line1P::Point(n);
				integral += point.W * f(point.u);
			}
		}
		std::cout << "integral = " << integral << std::endl;
	}
}


// Точка Гаусса - структура. Схема интегрирования - безымянная структура.
// При смене числа точек требует правки в двух местах: заголовок цикла и строка суммирования
namespace Test12 {
	namespace Gauss {
		struct Point1D {
			double u;
			double W;
		};

		const struct {
			int Points;

			Point1D Point(int n) const
			{
				static Point1D points[] = {0.0, 2.0};
				return points[n];
			}
		} Line1P = {1};

	}

	void run() {
		double integral = 0;
		for(int i = 0; i < N; i++) {
			for(int n = 0; n < Gauss::Line1P.Points; n++) {
				Gauss::Point1D point = Gauss::Line1P.Point(n);
				integral += point.W * f(point.u);
			}
		}
		std::cout << "integral = " << integral << std::endl;
	}
}

// Точка Гаусса - структура. Схема интегрирования - структура.
// При смене числа точек требует единственного исправления - строки получения схемы
// Пространство имен замусорено посторонними функциями
// Тормозное решение
namespace Test13 {
	namespace Gauss {
		struct Point1D {
			double u;
			double W;
		};

		struct Scheme {
			int Points;
			Point1D (*Point)(int n);
		};

		Point1D Line1P_Point(int n)
		{
			static Point1D points[] = {0.0, 2.0};
			return points[n];
		}

		const Scheme Line1P = {1, Line1P_Point};
	}

	void run() {
		double integral = 0;
		for(int i = 0; i < N; i++) {
			Gauss::Scheme scheme = Gauss::Line1P;
			for(int n = 0; n < scheme.Points; n++) {
				Gauss::Point1D point = scheme.Point(n);
				integral += point.W * f(point.u);
			}
		}
		std::cout << "integral = " << integral << std::endl;
	}
}


// Точка Гаусса - класс. Схема интегрирования - класс.
// При смене числа точек требует единственного исправления - строки получения схемы
namespace Test14 {
	namespace Gauss {
		class Point1D {
		public:
			Point1D(double u, double Weight): m_u(u), m_Weight(Weight) {}
			double u() { return m_u; }
			double Weight() { return m_Weight; }
		private:
			double m_u;
			double m_Weight;
		};

		namespace Scheme {
			Point1D Line_1_Points[] = {
				Point1D(0.0, 2.0)
			};

			class Line1 {
			public:
				int Points() {return 1; }
				Point1D Point(int n) { return Line_1_Points[n]; }
			};

		}

	}

	void run() {
		double integral = 0;
		for(int i = 0; i < N; i++) {
			Gauss::Scheme::Line1 scheme;
			for(int n = 0; n < scheme.Points(); n++) {
				Gauss::Point1D point = scheme.Point(n);
				integral += point.Weight() * f(point.u());
			}
		}
		std::cout << "integral = " << integral << std::endl;
	}
}


// Точка Гаусса - класс. Схема интегрирования - универсальный класс с массивом точек Гаусса.
// При смене числа точек требует единственного исправления - строки получения схемы
namespace Test15 {
	namespace Gauss {
		class Point1D {
		public:
			Point1D(double u, double Weight): m_u(u), m_Weight(Weight) {}
			double u() { return m_u; }
			double Weight() { return m_Weight; }
		private:
			double m_u;
			double m_Weight;
		};

		namespace Scheme {
			Point1D Line_1_Points[] = {
				Point1D(0.0, 2.0)
			};

			Point1D* Line_Points[] = {
				&Line_1_Points[0]
			};

			class Line {
			public:
				Line(int nPoints): nPoints(nPoints), points(Line_Points[nPoints-1]) {}
				int Points() {return nPoints; }
				Point1D Point(int n) { return points[n]; }
			private:
				int nPoints;
				Point1D* points;
			};

		}

	}

	void run() {
		double integral = 0;
		for(int i = 0; i < N; i++) {
			Gauss::Scheme::Line scheme(1);
			for(int n = 0; n < scheme.Points(); n++) {
				Gauss::Point1D point = scheme.Point(n);
				integral += point.Weight() * f(point.u());
			}
		}
		std::cout << "integral = " << integral << std::endl;
	}
}


int main()
{
	Timer::Start("Test 1");
	Test1::run();
	Timer::Stop();

	Timer::Start("Test 2");
	Test2::run();
	Timer::Stop();

	Timer::Start("Test 3");
	Test3::run();
	Timer::Stop();

	Timer::Start("Test 4");
	Test4::run();
	Timer::Stop();

	Timer::Start("Test 5");
	Test5::run();
	Timer::Stop();

	Timer::Start("Test 6");
	Test6::run();
	Timer::Stop();

	Timer::Start("Test 7");
	Test7::run();
	Timer::Stop();

	Timer::Start("Test 8");
	Test8::run();
	Timer::Stop();

	Timer::Start("Test 9");
	Test9::run();
	Timer::Stop();

	Timer::Start("Test 10");
	Test10::run();
	Timer::Stop();

	Timer::Start("Test 11");
	Test11::run();
	Timer::Stop();

	Timer::Start("Test 12");
	Test12::run();
	Timer::Stop();

	Timer::Start("Test 13");
	Test13::run();
	Timer::Stop();

	Timer::Start("Test 14");
	Test14::run();
	Timer::Stop();

	Timer::Start("Test 15");
	Test15::run();
	Timer::Stop();
}

Я удивлен, что стандартные решения через указатели на функции (стиль C) и через наследование (стиль С++) - два единственных тормозных способа упаковать константы в какую-то легко используемую оболочку; а городильня с массивами и универсальными классами равна по скорости прямому указанию констант в коде. Даже интересно, сколько всего мы еще по незнанию наделаем не так как надо .

Кстати, есть вопрос. В коде GetDP для численного интегрирования прошиты координаты точек Гаусса до 20-точечных схем включительно. Кто знает, откуда они взяли формулы? В литературе только пятиточечные схемы нахожу. В коде GetDP в последних знаках погрешности видны, в стиле "5.55555555555555552" - получается, как в анекдоте, неаккуратненько .

[zamtmn]

Если записывать данные не "влоб" - разница будет минимальной.
Например прокешировав все\частично в озу, а не пологаясь полностью на ос

[Нубий-IV]

Сообщение от zamtmn прокешировав все\частично в озу, а не пологаясь полностью на ос

В тесте стандартный кешированный поток C++ на порядок проиграл по скорости системной функции записи в файл. Еще Гугель утверждает, что отображение файлов в память даже промежуточный буфер не использует; а если памяти хватает, то и собственно файл не записывается - только однократное чтение/запись в память происходит: пишет передающая программа, читает принимающая. Как понимаю, это стандартное системное кеширование из коробки, желающие могут взамен написать свое. Опять же, не разбираюсь в тонкостях, все надо тестировать.

[zamtmn]

незнаю насчет c++, на паскале куча вариантов реализаций. сомневаюсь что в них можно будет найти проигрывающие стандартно-минимальному способу

Если начинать из такого долека, результат будет лет через 15-20))

[Нубий-IV]

Оказывается, гаусс-то засекречен! В книгах по интегрированию либо теория с примерами для низких степеней, либо числовые коэффициенты с 5-15 знаками, что недотягивает до точности double. Изредка попадаются честные двадцатизначные таблицы или формулы.
Книги:

• Приближенное вычисление интегралов.1967.Крылов
• Интерполяционные кубатурные формулы.1981.Мысовских
• Approximate calculation of multiple integrals.1971.Stroud
• Геометрическое моделирование.2002.Голованов

Что нашел - набрал. Для отрезков - 9 степень, для треугольников - 5, для прямоугольников - 9. Если понадобятся степени выше, можно уточнить книжные коэффициенты в максиме, численно решив нелинейную систему уравнений с начальным приближением.

Код:

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

 
#ifndef GAUSS_H
#define GAUSS_H


namespace Gauss {


class Point1D {
public:
	Point1D(double u, double Weight);
	double u()      const;
	double Weight() const;
private:
	double m_u;
	double m_Weight;
};


class Point2D {
public:
	Point2D(double u, double v, double Weight);
	double u()      const;
	double v()      const;
	double Weight() const;
private:
	double m_u;
	double m_v;
	double m_Weight;
};


class Point3D {
public:
	Point3D(double u, double v, double w, double Weight);
	double u()      const;
	double v()      const;
	double w()      const;
	double Weight() const;
private:
	double m_u;
	double m_v;
	double m_w;
	double m_Weight;
};


namespace Scheme {


class Line_1D_1P        { public: int NumPoints() const; Point1D Point(int n) const; private: static const Point1D points[]; };
class Line_3D_2P        { public: int NumPoints() const; Point1D Point(int n) const; private: static const Point1D points[]; };
class Line_5D_3P        { public: int NumPoints() const; Point1D Point(int n) const; private: static const Point1D points[]; };
class Line_7D_4P        { public: int NumPoints() const; Point1D Point(int n) const; private: static const Point1D points[]; };
class Line_9D_5P        { public: int NumPoints() const; Point1D Point(int n) const; private: static const Point1D points[]; };


class Triangle_1D_1P    { public: int NumPoints() const; Point2D Point(int n) const; private: static const Point2D points[]; };
class Triangle_2D_3P    { public: int NumPoints() const; Point2D Point(int n) const; private: static const Point2D points[]; };
class Triangle_3D_4P    { public: int NumPoints() const; Point2D Point(int n) const; private: static const Point2D points[]; };
class Triangle_4D_6P    { public: int NumPoints() const; Point2D Point(int n) const; private: static const Point2D points[]; };
class Triangle_5D_7P    { public: int NumPoints() const; Point2D Point(int n) const; private: static const Point2D points[]; };


class Quadrangle_1D_1P  { public: int NumPoints() const; Point2D Point(int n) const; private: static const Point2D points[]; };
class Quadrangle_2D_3P  { public: int NumPoints() const; Point2D Point(int n) const; private: static const Point2D points[]; };
class Quadrangle_3D_4P  { public: int NumPoints() const; Point2D Point(int n) const; private: static const Point2D points[]; };
class Quadrangle_5D_7P  { public: int NumPoints() const; Point2D Point(int n) const; private: static const Point2D points[]; };
class Quadrangle_7D_12P { public: int NumPoints() const; Point2D Point(int n) const; private: static const Point2D points[]; };
class Quadrangle_9D_17P { public: int NumPoints() const; Point2D Point(int n) const; private: static const Point2D points[]; };


} /* namespace Scheme */


} /* namespace Gauss */


#endif /* GAUSS_H */

Код:

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

 
#include "Gauss.h"


namespace Gauss {


Point1D::Point1D(double u, double Weight): m_u(u), m_Weight(Weight) {}
double Point1D::u()      const { return m_u; }
double Point1D::Weight() const { return m_Weight; }


Point2D::Point2D(double u, double v, double Weight): m_u(u), m_v(v), m_Weight(Weight) {}
double Point2D::u()      const { return m_u; }
double Point2D::v()      const { return m_v; }
double Point2D::Weight() const { return m_Weight; }


Point3D::Point3D(double u, double v, double w, double Weight): m_u(u), m_v(v), m_w(w), m_Weight(Weight) {}
double Point3D::u()      const { return m_u; }
double Point3D::v()      const { return m_v; }
double Point3D::w()      const { return m_w; }
double Point3D::Weight() const { return m_Weight; }


namespace Scheme {


int Line_1D_1P::NumPoints() const { return 1; }
int Line_3D_2P::NumPoints() const { return 2; }
int Line_5D_3P::NumPoints() const { return 3; }
int Line_7D_4P::NumPoints() const { return 4; }
int Line_9D_5P::NumPoints() const { return 5; }


Point1D Line_1D_1P::Point(int n) const { return points[n]; }
Point1D Line_3D_2P::Point(int n) const { return points[n]; }
Point1D Line_5D_3P::Point(int n) const { return points[n]; }
Point1D Line_7D_4P::Point(int n) const { return points[n]; }
Point1D Line_9D_5P::Point(int n) const { return points[n]; }


const Point1D Line_1D_1P::points[] = {
	Point1D( 0.0000000000000000e+00, 2.0000000000000000e+00),
};

const Point1D Line_3D_2P::points[] = {
	Point1D(-5.7735026918962576e-01, 1.0000000000000000e+00),
	Point1D( 5.7735026918962576e-01, 1.0000000000000000e+00),
};

const Point1D Line_5D_3P::points[] = {
	Point1D(-7.7459666924148338e-01, 5.5555555555555556e-01),
	Point1D( 0.0000000000000000e+00, 8.8888888888888889e-01),
	Point1D( 7.7459666924148338e-01, 5.5555555555555556e-01),
};

const Point1D Line_7D_4P::points[] = {
	Point1D(-8.6113631159405258e-01, 3.4785484513745386e-01),
	Point1D(-3.3998104358485626e-01, 6.5214515486254614e-01),
	Point1D( 3.3998104358485626e-01, 6.5214515486254614e-01),
	Point1D( 8.6113631159405258e-01, 3.4785484513745386e-01),
};

const Point1D Line_9D_5P::points[] = {
	Point1D(-9.0617984593866399e-01, 2.3692688505618909e-01),
	Point1D(-5.3846931010568309e-01, 4.7862867049936647e-01),
	Point1D( 0.0000000000000000e+00, 5.6888888888888889e-01),
	Point1D( 5.3846931010568309e-01, 4.7862867049936647e-01),
	Point1D( 9.0617984593866399e-01, 2.3692688505618909e-01),
};


int Triangle_1D_1P::NumPoints() const { return 1; }
int Triangle_2D_3P::NumPoints() const { return 3; }
int Triangle_3D_4P::NumPoints() const { return 4; }
int Triangle_4D_6P::NumPoints() const { return 6; }
int Triangle_5D_7P::NumPoints() const { return 7; }


Point2D Triangle_1D_1P::Point(int n) const { return points[n]; }
Point2D Triangle_2D_3P::Point(int n) const { return points[n]; }
Point2D Triangle_3D_4P::Point(int n) const { return points[n]; }
Point2D Triangle_4D_6P::Point(int n) const { return points[n]; }
Point2D Triangle_5D_7P::Point(int n) const { return points[n]; }


const Point2D Triangle_1D_1P::points[] = {
	Point2D( 3.3333333333333333e-01,  3.3333333333333333e-01, 5.0000000000000000e-01),
};

const Point2D Triangle_2D_3P::points[] = {
	Point2D( 1.6666666666666667e-01,  6.6666666666666667e-01, 1.6666666666666667e-01),
	Point2D( 1.6666666666666667e-01,  1.6666666666666667e-01, 1.6666666666666667e-01),
	Point2D( 6.6666666666666667e-01,  1.6666666666666667e-01, 1.6666666666666667e-01),
};

const Point2D Triangle_3D_4P::points[] = {
	Point2D( 2.0000000000000000e-01,  2.0000000000000000e-01, 2.6041666666666667e-01),
	Point2D( 3.3333333333333333e-01,  3.3333333333333333e-01,-2.8125000000000000e-01),
	Point2D( 2.0000000000000000e-01,  6.0000000000000000e-01, 2.6041666666666667e-01),
	Point2D( 6.0000000000000000e-01,  2.0000000000000000e-01, 2.6041666666666667e-01),
};

const Point2D Triangle_4D_6P::points[] = {
	Point2D( 8.1684757298045851e-01, 9.1576213509770743e-02, 5.4975871827660934e-02),
	Point2D( 9.1576213509770743e-02, 8.1684757298045851e-01, 5.4975871827660934e-02),
	Point2D( 9.1576213509770743e-02, 9.1576213509770743e-02, 5.4975871827660934e-02),
	Point2D( 1.0810301816807023e-01, 4.4594849091596489e-01, 1.1169079483900573e-01),
	Point2D( 4.4594849091596489e-01, 1.0810301816807023e-01, 1.1169079483900573e-01),
	Point2D( 4.4594849091596489e-01, 4.4594849091596489e-01, 1.1169079483900573e-01),
};

const Point2D Triangle_5D_7P::points[] = {
	Point2D( 1.0128650732345634e-01,  1.0128650732345634e-01, 6.2969590272413576e-02),
	Point2D( 1.0128650732345634e-01,  7.9742698535308732e-01, 6.2969590272413576e-02),
	Point2D( 3.3333333333333333e-01,  3.3333333333333333e-01, 1.1250000000000000e-01),
	Point2D( 4.7014206410511509e-01,  4.7014206410511509e-01, 6.6197076394253090e-02),
	Point2D( 4.7014206410511509e-01,  5.9715871789769820e-02, 6.6197076394253090e-02),
	Point2D( 5.9715871789769820e-02,  4.7014206410511509e-01, 6.6197076394253090e-02),
	Point2D( 7.9742698535308732e-01,  1.0128650732345634e-01, 6.2969590272413576e-02),
};


int Quadrangle_1D_1P ::NumPoints() const { return  1; }
int Quadrangle_2D_3P ::NumPoints() const { return  3; }
int Quadrangle_3D_4P ::NumPoints() const { return  4; }
int Quadrangle_5D_7P ::NumPoints() const { return  7; }
int Quadrangle_7D_12P::NumPoints() const { return 12; }
int Quadrangle_9D_17P::NumPoints() const { return 17; }


Point2D Quadrangle_1D_1P ::Point(int n) const { return points[n]; }
Point2D Quadrangle_2D_3P ::Point(int n) const { return points[n]; }
Point2D Quadrangle_3D_4P ::Point(int n) const { return points[n]; }
Point2D Quadrangle_5D_7P ::Point(int n) const { return points[n]; }
Point2D Quadrangle_7D_12P::Point(int n) const { return points[n]; }
Point2D Quadrangle_9D_17P::Point(int n) const { return points[n]; }


const Point2D Quadrangle_1D_1P::points[] = {
	Point2D( 0.0000000000000000e+00,  0.0000000000000000e+00, 4.0000000000000000e+00),
};

const Point2D Quadrangle_2D_3P::points[] = {
	Point2D(-4.0824829046386302e-01,-7.0710678118654752e-01, 1.3333333333333333e+00),
	Point2D(-4.0824829046386302e-01, 7.0710678118654752e-01, 1.3333333333333333e+00),
	Point2D( 8.1649658092772603e-01, 0.0000000000000000e+00, 1.3333333333333333e+00),

};

const Point2D Quadrangle_3D_4P::points[] = {
	Point2D(-5.7735026918962576e-01, -5.7735026918962576e-01, 1.0000000000000000e+00),
	Point2D(-5.7735026918962576e-01,  5.7735026918962576e-01, 1.0000000000000000e+00),
	Point2D( 5.7735026918962576e-01, -5.7735026918962576e-01, 1.0000000000000000e+00),
	Point2D( 5.7735026918962576e-01,  5.7735026918962576e-01, 1.0000000000000000e+00),
};

const Point2D Quadrangle_5D_7P::points[] = {
	Point2D(-7.7459666924148338e-01, -5.7735026918962576e-01, 5.5555555555555556e-01),
	Point2D(-7.7459666924148338e-01,  5.7735026918962576e-01, 5.5555555555555556e-01),
	Point2D( 0.0000000000000000e+00, -9.6609178307929590e-01, 3.1746031746031746e-01),
	Point2D( 0.0000000000000000e+00,  0.0000000000000000e+00, 1.1428571428571429e+00),
	Point2D( 0.0000000000000000e+00,  9.6609178307929590e-01, 3.1746031746031746e-01),
	Point2D( 7.7459666924148338e-01, -5.7735026918962576e-01, 5.5555555555555556e-01),
	Point2D( 7.7459666924148338e-01,  5.7735026918962576e-01, 5.5555555555555556e-01),
};

const Point2D Quadrangle_7D_12P::points[] = {
	Point2D(-9.2582009977255146e-01,  0.0000000000000000e+00, 2.4197530864197531e-01),
	Point2D(-8.0597978291859874e-01, -8.0597978291859874e-01, 2.3743177469063023e-01),
	Point2D(-8.0597978291859874e-01,  8.0597978291859874e-01, 2.3743177469063023e-01),
	Point2D(-3.8055443320831566e-01, -3.8055443320831566e-01, 5.2059291666739446e-01),
	Point2D(-3.8055443320831566e-01,  3.8055443320831566e-01, 5.2059291666739446e-01),
	Point2D( 0.0000000000000000e+00, -9.2582009977255146e-01, 2.4197530864197531e-01),
	Point2D( 0.0000000000000000e+00,  9.2582009977255146e-01, 2.4197530864197531e-01),
	Point2D( 3.8055443320831566e-01, -3.8055443320831566e-01, 5.2059291666739446e-01),
	Point2D( 3.8055443320831566e-01,  3.8055443320831566e-01, 5.2059291666739446e-01),
	Point2D( 8.0597978291859874e-01, -8.0597978291859874e-01, 2.3743177469063023e-01),
	Point2D( 8.0597978291859874e-01,  8.0597978291859874e-01, 2.3743177469063023e-01),
	Point2D( 9.2582009977255146e-01,  0.0000000000000000e+00, 2.4197530864197531e-01),
};

const Point2D Quadrangle_9D_17P::points[] = {
	Point2D(-9.6884996636197772e-01, -6.3068011973166885e-01, 8.8879378170198707e-02),
	Point2D(-9.2796164595956967e-01,  7.5027709997890053e-01, 1.1209960212959649e-01),
	Point2D(-8.5261572933366231e-01,  7.6208328192617173e-02, 2.6905133763978080e-01),
	Point2D(-7.6208328192617173e-02, -8.5261572933366231e-01, 2.6905133763978080e-01),
	Point2D(-7.5027709997890053e-01, -9.2796164595956967e-01, 1.1209960212959649e-01),
	Point2D(-6.3068011973166885e-01,  9.6884996636197772e-01, 8.8879378170198707e-02),
	Point2D(-5.2373582021442934e-01, -4.5333982113564719e-01, 3.9828243926207010e-01),
	Point2D(-4.5333982113564719e-01,  5.2373582021442934e-01, 3.9828243926207010e-01),
	Point2D( 0.0000000000000000e+00,  0.0000000000000000e+00, 5.2674897119341564e-01),
	Point2D( 4.5333982113564719e-01, -5.2373582021442934e-01, 3.9828243926207010e-01),
	Point2D( 5.2373582021442934e-01,  4.5333982113564719e-01, 3.9828243926207010e-01),
	Point2D( 6.3068011973166885e-01, -9.6884996636197772e-01, 8.8879378170198707e-02),
	Point2D( 7.5027709997890053e-01,  9.2796164595956967e-01, 1.1209960212959649e-01),
	Point2D( 7.6208328192617173e-02,  8.5261572933366231e-01, 2.6905133763978080e-01),
	Point2D( 8.5261572933366231e-01, -7.6208328192617173e-02, 2.6905133763978080e-01),
	Point2D( 9.2796164595956967e-01, -7.5027709997890053e-01, 1.1209960212959649e-01),
	Point2D( 9.6884996636197772e-01,  6.3068011973166885e-01, 8.8879378170198707e-02),
};


} /* namespace Scheme */


} /* namespace Gauss */

Пример использования (Quadrangle_3D_4P значит "Прямоугольник", "Полином до 3 степени", "4 точки"):

Код:

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

 
#include<iostream>
#include "Gauss.h"


double f(double u, double v)
{
	return u*u + u*v + v*v;
}


int main()
{
	double Integral = 0.0;

	Gauss::Scheme::Quadrangle_3D_4P scheme;

	for(int n = 0; n < scheme.NumPoints(); n++) {
		Gauss::Point2D point = scheme.Point(n);
		Integral += point.Weight() * f(point.u(), point.v());
	}

	std::cout << "Integral = " << Integral << std::endl;

	return 0;
}

Тест точности:

Код:

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

 
#include<iostream>
#include<ctime>
#include<math.h>
#include "Gauss.h"


struct Line_Poly_1D {
	double f(double u) {
		return
		1 +
		u;
	}
	double F() { return 2.0; }
};

struct Line_Poly_2D {
	double f(double u) {
		return
		1 +
		u +
		u*u;
	}
	double F() { return 8.0 / 3.0; }
};

struct Line_Poly_3D {
	double f(double u){
		return
		1 +
		u +
		u*u +
		u*u*u;
	}
	double F() { return 8.0 / 3.0; }
};

struct Line_Poly_4D {
	double f(double u) {
		return
		1 +
		u +
		u*u +
		u*u*u +
		u*u*u*u;
	}
	double F() { return 46.0 / 15.0; }
};

struct Line_Poly_5D {
	double f(double u) {
		return
		1 +
		u +
		u*u +
		u*u*u +
		u*u*u*u +
		u*u*u*u*u;
	}
	double F() { return 46.0 / 15.0; }
};

struct Line_Poly_6D {
	double f(double u) {
		return
		1 +
		u +
		u*u +
		u*u*u +
		u*u*u*u +
		u*u*u*u*u +
		u*u*u*u*u*u;
	}
	double F() { return 352.0 / 105.0; }
};

struct Line_Poly_7D {
	double f(double u) {
		return
		1 +
		u +
		u*u +
		u*u*u +
		u*u*u*u +
		u*u*u*u*u +
		u*u*u*u*u*u +
		u*u*u*u*u*u*u;
	}
	double F() { return 352.0 / 105.0; }
};

struct Line_Poly_8D {
	double f(double u) {
		return
		1 +
		u +
		u*u +
		u*u*u +
		u*u*u*u +
		u*u*u*u*u +
		u*u*u*u*u*u +
		u*u*u*u*u*u*u +
		u*u*u*u*u*u*u*u;
	}
	double F() { return 1126.0 / 315.0; }
};

struct Line_Poly_9D {
	double f(double u) {
		return
		1 +
		u +
		u*u +
		u*u*u +
		u*u*u*u +
		u*u*u*u*u +
		u*u*u*u*u*u +
		u*u*u*u*u*u*u +
		u*u*u*u*u*u*u*u +
		u*u*u*u*u*u*u*u*u;
	}
	double F() { return 1126.0 / 315.0; }
};

#define TEST_LINE(POLYNOM, GAUSS_SCHEME) \
{ \
	std::cout.precision(17); \
	std::cout << "Gauss::Scheme::" << #GAUSS_SCHEME  << std::endl; \
	std::cout << #POLYNOM << std::endl; \
	POLYNOM poly; \
	Gauss::Scheme::GAUSS_SCHEME scheme; \
	double F = 0.0; \
	for(int n = 0; n < scheme.NumPoints(); n++) { \
		Gauss::Point1D point = scheme.Point(n); \
		F += point.Weight() * poly.f(point.u()); \
	} \
	std::cout << "Numerically:  "  <<      F   << std::endl; \
	std::cout << "Analytically: " << poly.F() << std::endl << std::endl; \
}


struct Triangle_Poly_1D {
	double f(double u, double v) {
		return
		1 +
		u + v;
	}
	double F() { return 5.0 / 6.0; }
};

struct Triangle_Poly_2D {
	double f(double u, double v) {
		return
		1 +
		u + v +
		u*u + v*v + u*v;
	}
	double F() { return 25.0 / 24.0; }
};

struct Triangle_Poly_3D {
	double f(double u, double v) {
		return
		1 +
		u + v +
		u*u + v*v + u*v +
		u*u*u + u*u*v + u*v*v + v*v*v;
	}
	double F() { return 47.0 / 40.0; }
};

struct Triangle_Poly_4D {
	double f(double u, double v) {
		return
		1 +
		u + v +
		u*u + v*v + u*v +
		u*u*u + u*u*v + u*v*v + v*v*v +
		u*u*u*u + u*u*u*v + u*u*v*v + u*v*v*v + v*v*v*v;
	}
	double F() { return 91.0 / 72.0; }
};

struct Triangle_Poly_5D {
	double f(double u, double v) {
		return
		1 + u + v +
		u*u + v*v + u*v +
		u*u*u + u*u*v + u*v*v + v*v*v +
		u*u*u*u + u*u*u*v + u*u*v*v + u*v*v*v + v*v*v*v +
		u*u*u*u*u + u*u*u*u*v + u*u*u*v*v + u*u*v*v*v + u*v*v*v*v + v*v*v*v*v;
	}
	double F() { return 3341.0 / 2520.0; }
};

#define TEST_TRIANGLE(POLYNOM, GAUSS_SCHEME) \
{ \
	std::cout.precision(17); \
	std::cout << "Gauss::Scheme::" << #GAUSS_SCHEME  << std::endl; \
	std::cout << #POLYNOM << std::endl; \
	POLYNOM poly; \
	Gauss::Scheme::GAUSS_SCHEME scheme; \
	double F = 0.0; \
	for(int n = 0; n < scheme.NumPoints(); n++) { \
		Gauss::Point2D point = scheme.Point(n); \
		F += point.Weight() * poly.f(point.u(), point.v()); \
	} \
	std::cout << "Numerically:  "  <<      F   << std::endl; \
	std::cout << "Analytically: " << poly.F() << std::endl << std::endl; \
}


struct Quadrangle_Poly_1D {
	double f(double u, double v) {
		return
		1 +
		u + v;
	}
	double F() { return 4.0; }
};

struct Quadrangle_Poly_2D {
	double f(double u, double v) {
		return
		1 +
		u + v +
		u*u + v*v + u*v;
	}
	double F() { return 20.0 / 3.0; }
};

struct Quadrangle_Poly_3D {
	double f(double u, double v) {
		return
		1 +
		u + v +
		u*u + v*v + u*v +
		u*u*u + u*u*v + u*v*v + v*v*v;
	}
	double F() { return 20.0 / 3.0; }
};

struct Quadrangle_Poly_4D {
	double f(double u, double v) {
		return
		1 +
		u + v +
		u*u + v*v + u*v +
		u*u*u + u*u*v + u*v*v + v*v*v +
		u*u*u*u + u*u*u*v + u*u*v*v + u*v*v*v + v*v*v*v;
	}
	double F() { return 392.0 / 45.0; }
};

struct Quadrangle_Poly_5D {
	double f(double u, double v) {
		return
		1 + u + v +
		u*u + v*v + u*v +
		u*u*u + u*u*v + u*v*v + v*v*v +
		u*u*u*u + u*u*u*v + u*u*v*v + u*v*v*v + v*v*v*v +
		u*u*u*u*u + u*u*u*u*v + u*u*u*v*v + u*u*v*v*v + u*v*v*v*v + v*v*v*v*v;
	}
	double F() { return 392.0 / 45.0; }
};

struct Quadrangle_Poly_6D {
	double f(double u, double v) {
		return
		1 + u + v +
		u*u + v*v + u*v +
		u*u*u + u*u*v + u*v*v + v*v*v +
		u*u*u*u + u*u*u*v + u*u*v*v + u*v*v*v + v*v*v*v +
		u*u*u*u*u + u*u*u*u*v + u*u*u*v*v + u*u*v*v*v + u*v*v*v*v + v*v*v*v*v +
		u*u*u*u*u*u + u*u*u*u*u*v + u*u*u*u*v*v + u*u*u*v*v*v + u*u*v*v*v*v + u*v*v*v*v*v + v*v*v*v*v*v;
	}
	double F() { return 3272.0 / 315.0; }
};


struct Quadrangle_Poly_7D {
	double f(double u, double v) {
		return
		1 + u + v +
		u*u + v*v + u*v +
		u*u*u + u*u*v + u*v*v + v*v*v +
		u*u*u*u + u*u*u*v + u*u*v*v + u*v*v*v + v*v*v*v +
		u*u*u*u*u + u*u*u*u*v + u*u*u*v*v + u*u*v*v*v + u*v*v*v*v + v*v*v*v*v +
		u*u*u*u*u*u + u*u*u*u*u*v + u*u*u*u*v*v + u*u*u*v*v*v + u*u*v*v*v*v + u*v*v*v*v*v + v*v*v*v*v*v +
		u*u*u*u*u*u*u + u*u*u*u*u*u*v + u*u*u*u*u*v*v + u*u*u*u*v*v*v + u*u*u*v*v*v*v + u*u*v*v*v*v*v + u*v*v*v*v*v*v + v*v*v*v*v*v*v;
	}
	double F() { return 3272.0 / 315.0; }
};

struct Quadrangle_Poly_8D {
	double f(double u, double v) {
		return
		1 + u + v +
		u*u + v*v + u*v +
		u*u*u + u*u*v + u*v*v + v*v*v +
		u*u*u*u + u*u*u*v + u*u*v*v + u*v*v*v + v*v*v*v +
		u*u*u*u*u + u*u*u*u*v + u*u*u*v*v + u*u*v*v*v + u*v*v*v*v + v*v*v*v*v +
		u*u*u*u*u*u + u*u*u*u*u*v + u*u*u*u*v*v + u*u*u*v*v*v + u*u*v*v*v*v + u*v*v*v*v*v + v*v*v*v*v*v +
		u*u*u*u*u*u*u + u*u*u*u*u*u*v + u*u*u*u*u*v*v + u*u*u*u*v*v*v + u*u*u*v*v*v*v + u*u*v*v*v*v*v + u*v*v*v*v*v*v + v*v*v*v*v*v*v +
		u*u*u*u*u*u*u*u + u*u*u*u*u*u*u*v + u*u*u*u*u*u*v*v + u*u*u*u*u*v*v*v + u*u*u*u*v*v*v*v + u*u*u*v*v*v*v*v + u*u*v*v*v*v*v*v + u*v*v*v*v*v*v*v + v*v*v*v*v*v*v*v;
	}
	double F() { return 2068.0 / 175.0; }
};

struct Quadrangle_Poly_9D {
	double f(double u, double v) {
		return
		1 + u + v +
		u*u + v*v + u*v +
		u*u*u + u*u*v + u*v*v + v*v*v +
		u*u*u*u + u*u*u*v + u*u*v*v + u*v*v*v + v*v*v*v +
		u*u*u*u*u + u*u*u*u*v + u*u*u*v*v + u*u*v*v*v + u*v*v*v*v + v*v*v*v*v +
		u*u*u*u*u*u + u*u*u*u*u*v + u*u*u*u*v*v + u*u*u*v*v*v + u*u*v*v*v*v + u*v*v*v*v*v + v*v*v*v*v*v +
		u*u*u*u*u*u*u + u*u*u*u*u*u*v + u*u*u*u*u*v*v + u*u*u*u*v*v*v + u*u*u*v*v*v*v + u*u*v*v*v*v*v + u*v*v*v*v*v*v + v*v*v*v*v*v*v +
		u*u*u*u*u*u*u*u + u*u*u*u*u*u*u*v + u*u*u*u*u*u*v*v + u*u*u*u*u*v*v*v + u*u*u*u*v*v*v*v + u*u*u*v*v*v*v*v + u*u*v*v*v*v*v*v + u*v*v*v*v*v*v*v + v*v*v*v*v*v*v*v +
		u*u*u*u*u*u*u*u*u + u*u*u*u*u*u*u*u*v + u*u*u*u*u*u*u*v*v + u*u*u*u*u*u*v*v*v + u*u*u*u*u*v*v*v*v + u*u*u*u*v*v*v*v*v + u*u*u*v*v*v*v*v*v + u*u*v*v*v*v*v*v*v + u*v*v*v*v*v*v*v*v + v*v*v*v*v*v*v*v*v;
	}
	double F() { return 2068.0 / 175.0; }
};

#define TEST_QUADRANGLE(POLYNOM, GAUSS_SCHEME) \
{ \
	std::cout.precision(17); \
	std::cout << "Gauss::Scheme::" << #GAUSS_SCHEME  << std::endl; \
	std::cout << #POLYNOM << std::endl; \
	POLYNOM poly; \
	Gauss::Scheme::GAUSS_SCHEME scheme; \
	double F = 0.0; \
	for(int n = 0; n < scheme.NumPoints(); n++) { \
		Gauss::Point2D point = scheme.Point(n); \
		F += point.Weight() * poly.f(point.u(), point.v()); \
	} \
	std::cout << "Numerically:  "  <<      F   << std::endl; \
	std::cout << "Analytically: " << poly.F() << std::endl << std::endl; \
}


int main()
{
	TEST_LINE(Line_Poly_1D, Line_1D_1P)

	TEST_LINE(Line_Poly_1D, Line_3D_2P)
	TEST_LINE(Line_Poly_2D, Line_3D_2P)
	TEST_LINE(Line_Poly_3D, Line_3D_2P)

	TEST_LINE(Line_Poly_1D, Line_5D_3P)
	TEST_LINE(Line_Poly_2D, Line_5D_3P)
	TEST_LINE(Line_Poly_3D, Line_5D_3P)
	TEST_LINE(Line_Poly_4D, Line_5D_3P)
	TEST_LINE(Line_Poly_5D, Line_5D_3P)

	TEST_LINE(Line_Poly_1D, Line_7D_4P)
	TEST_LINE(Line_Poly_2D, Line_7D_4P)
	TEST_LINE(Line_Poly_3D, Line_7D_4P)
	TEST_LINE(Line_Poly_4D, Line_7D_4P)
	TEST_LINE(Line_Poly_5D, Line_7D_4P)
	TEST_LINE(Line_Poly_6D, Line_7D_4P)
	TEST_LINE(Line_Poly_7D, Line_7D_4P)

	TEST_LINE(Line_Poly_1D, Line_9D_5P)
	TEST_LINE(Line_Poly_2D, Line_9D_5P)
	TEST_LINE(Line_Poly_3D, Line_9D_5P)
	TEST_LINE(Line_Poly_4D, Line_9D_5P)
	TEST_LINE(Line_Poly_5D, Line_9D_5P)
	TEST_LINE(Line_Poly_6D, Line_9D_5P)
	TEST_LINE(Line_Poly_7D, Line_9D_5P)
	TEST_LINE(Line_Poly_8D, Line_9D_5P)
	TEST_LINE(Line_Poly_9D, Line_9D_5P)


	TEST_TRIANGLE(Triangle_Poly_1D, Triangle_1D_1P)

	TEST_TRIANGLE(Triangle_Poly_1D, Triangle_2D_3P)
	TEST_TRIANGLE(Triangle_Poly_2D, Triangle_2D_3P)

	TEST_TRIANGLE(Triangle_Poly_1D, Triangle_3D_4P)
	TEST_TRIANGLE(Triangle_Poly_2D, Triangle_3D_4P)
	TEST_TRIANGLE(Triangle_Poly_3D, Triangle_3D_4P)

	TEST_TRIANGLE(Triangle_Poly_1D, Triangle_4D_6P)
	TEST_TRIANGLE(Triangle_Poly_2D, Triangle_4D_6P)
	TEST_TRIANGLE(Triangle_Poly_3D, Triangle_4D_6P)
	TEST_TRIANGLE(Triangle_Poly_4D, Triangle_4D_6P)

	TEST_TRIANGLE(Triangle_Poly_1D, Triangle_5D_7P)
	TEST_TRIANGLE(Triangle_Poly_2D, Triangle_5D_7P)
	TEST_TRIANGLE(Triangle_Poly_3D, Triangle_5D_7P)
	TEST_TRIANGLE(Triangle_Poly_4D, Triangle_5D_7P)
	TEST_TRIANGLE(Triangle_Poly_5D, Triangle_5D_7P)


	TEST_QUADRANGLE(Quadrangle_Poly_1D, Quadrangle_1D_1P)

	TEST_QUADRANGLE(Quadrangle_Poly_1D, Quadrangle_2D_3P)
	TEST_QUADRANGLE(Quadrangle_Poly_2D, Quadrangle_2D_3P)

	TEST_QUADRANGLE(Quadrangle_Poly_1D, Quadrangle_3D_4P)
	TEST_QUADRANGLE(Quadrangle_Poly_2D, Quadrangle_3D_4P)
	TEST_QUADRANGLE(Quadrangle_Poly_3D, Quadrangle_3D_4P)

	TEST_QUADRANGLE(Quadrangle_Poly_1D, Quadrangle_5D_7P)
	TEST_QUADRANGLE(Quadrangle_Poly_2D, Quadrangle_5D_7P)
	TEST_QUADRANGLE(Quadrangle_Poly_3D, Quadrangle_5D_7P)
	TEST_QUADRANGLE(Quadrangle_Poly_4D, Quadrangle_5D_7P)
	TEST_QUADRANGLE(Quadrangle_Poly_5D, Quadrangle_5D_7P)

	TEST_QUADRANGLE(Quadrangle_Poly_1D, Quadrangle_7D_12P)
	TEST_QUADRANGLE(Quadrangle_Poly_2D, Quadrangle_7D_12P)
	TEST_QUADRANGLE(Quadrangle_Poly_3D, Quadrangle_7D_12P)
	TEST_QUADRANGLE(Quadrangle_Poly_4D, Quadrangle_7D_12P)
	TEST_QUADRANGLE(Quadrangle_Poly_5D, Quadrangle_7D_12P)
	TEST_QUADRANGLE(Quadrangle_Poly_6D, Quadrangle_7D_12P)
	TEST_QUADRANGLE(Quadrangle_Poly_7D, Quadrangle_7D_12P)

	TEST_QUADRANGLE(Quadrangle_Poly_1D, Quadrangle_9D_17P)
	TEST_QUADRANGLE(Quadrangle_Poly_2D, Quadrangle_9D_17P)
	TEST_QUADRANGLE(Quadrangle_Poly_3D, Quadrangle_9D_17P)
	TEST_QUADRANGLE(Quadrangle_Poly_4D, Quadrangle_9D_17P)
	TEST_QUADRANGLE(Quadrangle_Poly_5D, Quadrangle_9D_17P)
	TEST_QUADRANGLE(Quadrangle_Poly_6D, Quadrangle_9D_17P)
	TEST_QUADRANGLE(Quadrangle_Poly_7D, Quadrangle_9D_17P)
	TEST_QUADRANGLE(Quadrangle_Poly_8D, Quadrangle_9D_17P)
	TEST_QUADRANGLE(Quadrangle_Poly_9D, Quadrangle_9D_17P)

	return 0;
}

Код:

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

Gauss::Scheme::Line_1D_1P
Line_Poly_1D
Numerically:  2
Analytically: 2

Gauss::Scheme::Line_3D_2P
Line_Poly_1D
Numerically:  2
Analytically: 2

Gauss::Scheme::Line_3D_2P
Line_Poly_2D
Numerically:  2.6666666666666665
Analytically: 2.6666666666666665

Gauss::Scheme::Line_3D_2P
Line_Poly_3D
Numerically:  2.6666666666666665
Analytically: 2.6666666666666665

Gauss::Scheme::Line_5D_3P
Line_Poly_1D
Numerically:  2
Analytically: 2

Gauss::Scheme::Line_5D_3P
Line_Poly_2D
Numerically:  2.666666666666667
Analytically: 2.6666666666666665

Gauss::Scheme::Line_5D_3P
Line_Poly_3D
Numerically:  2.666666666666667
Analytically: 2.6666666666666665

Gauss::Scheme::Line_5D_3P
Line_Poly_4D
Numerically:  3.0666666666666669
Analytically: 3.0666666666666669

Gauss::Scheme::Line_5D_3P
Line_Poly_5D
Numerically:  3.0666666666666669
Analytically: 3.0666666666666669

Gauss::Scheme::Line_7D_4P
Line_Poly_1D
Numerically:  1.9999999999999998
Analytically: 2

Gauss::Scheme::Line_7D_4P
Line_Poly_2D
Numerically:  2.6666666666666665
Analytically: 2.6666666666666665

Gauss::Scheme::Line_7D_4P
Line_Poly_3D
Numerically:  2.6666666666666665
Analytically: 2.6666666666666665

Gauss::Scheme::Line_7D_4P
Line_Poly_4D
Numerically:  3.0666666666666664
Analytically: 3.0666666666666669

Gauss::Scheme::Line_7D_4P
Line_Poly_5D
Numerically:  3.0666666666666664
Analytically: 3.0666666666666669

Gauss::Scheme::Line_7D_4P
Line_Poly_6D
Numerically:  3.3523809523809525
Analytically: 3.3523809523809525

Gauss::Scheme::Line_7D_4P
Line_Poly_7D
Numerically:  3.3523809523809529
Analytically: 3.3523809523809525

Gauss::Scheme::Line_9D_5P
Line_Poly_1D
Numerically:  2
Analytically: 2

Gauss::Scheme::Line_9D_5P
Line_Poly_2D
Numerically:  2.666666666666667
Analytically: 2.6666666666666665

Gauss::Scheme::Line_9D_5P
Line_Poly_3D
Numerically:  2.6666666666666665
Analytically: 2.6666666666666665

Gauss::Scheme::Line_9D_5P
Line_Poly_4D
Numerically:  3.0666666666666664
Analytically: 3.0666666666666669

Gauss::Scheme::Line_9D_5P
Line_Poly_5D
Numerically:  3.0666666666666664
Analytically: 3.0666666666666669

Gauss::Scheme::Line_9D_5P
Line_Poly_6D
Numerically:  3.352380952380952
Analytically: 3.3523809523809525

Gauss::Scheme::Line_9D_5P
Line_Poly_7D
Numerically:  3.352380952380952
Analytically: 3.3523809523809525

Gauss::Scheme::Line_9D_5P
Line_Poly_8D
Numerically:  3.5746031746031743
Analytically: 3.5746031746031748

Gauss::Scheme::Line_9D_5P
Line_Poly_9D
Numerically:  3.5746031746031743
Analytically: 3.5746031746031748

Gauss::Scheme::Triangle_1D_1P
Triangle_Poly_1D
Numerically:  0.83333333333333326
Analytically: 0.83333333333333337

Gauss::Scheme::Triangle_2D_3P
Triangle_Poly_1D
Numerically:  0.83333333333333326
Analytically: 0.83333333333333337

Gauss::Scheme::Triangle_2D_3P
Triangle_Poly_2D
Numerically:  1.0416666666666665
Analytically: 1.0416666666666667

Gauss::Scheme::Triangle_3D_4P
Triangle_Poly_1D
Numerically:  0.83333333333333348
Analytically: 0.83333333333333337

Gauss::Scheme::Triangle_3D_4P
Triangle_Poly_2D
Numerically:  1.0416666666666667
Analytically: 1.0416666666666667

Gauss::Scheme::Triangle_3D_4P
Triangle_Poly_3D
Numerically:  1.1750000000000003
Analytically: 1.175

Gauss::Scheme::Triangle_4D_6P
Triangle_Poly_1D
Numerically:  0.83333333333333326
Analytically: 0.83333333333333337

Gauss::Scheme::Triangle_4D_6P
Triangle_Poly_2D
Numerically:  1.0416666666666667
Analytically: 1.0416666666666667

Gauss::Scheme::Triangle_4D_6P
Triangle_Poly_3D
Numerically:  1.175
Analytically: 1.175

Gauss::Scheme::Triangle_4D_6P
Triangle_Poly_4D
Numerically:  1.2638888888888888
Analytically: 1.2638888888888888

Gauss::Scheme::Triangle_5D_7P
Triangle_Poly_1D
Numerically:  0.83333333333333326
Analytically: 0.83333333333333337

Gauss::Scheme::Triangle_5D_7P
Triangle_Poly_2D
Numerically:  1.0416666666666665
Analytically: 1.0416666666666667

Gauss::Scheme::Triangle_5D_7P
Triangle_Poly_3D
Numerically:  1.175
Analytically: 1.175

Gauss::Scheme::Triangle_5D_7P
Triangle_Poly_4D
Numerically:  1.2638888888888891
Analytically: 1.2638888888888888

Gauss::Scheme::Triangle_5D_7P
Triangle_Poly_5D
Numerically:  1.325793650793651
Analytically: 1.3257936507936507

Gauss::Scheme::Quadrangle_1D_1P
Quadrangle_Poly_1D
Numerically:  4
Analytically: 4

Gauss::Scheme::Quadrangle_2D_3P
Quadrangle_Poly_1D
Numerically:  3.9999999999999996
Analytically: 4

Gauss::Scheme::Quadrangle_2D_3P
Quadrangle_Poly_2D
Numerically:  6.666666666666667
Analytically: 6.666666666666667

Gauss::Scheme::Quadrangle_3D_4P
Quadrangle_Poly_1D
Numerically:  4
Analytically: 4

Gauss::Scheme::Quadrangle_3D_4P
Quadrangle_Poly_2D
Numerically:  6.6666666666666661
Analytically: 6.666666666666667

Gauss::Scheme::Quadrangle_3D_4P
Quadrangle_Poly_3D
Numerically:  6.6666666666666661
Analytically: 6.666666666666667

Gauss::Scheme::Quadrangle_5D_7P
Quadrangle_Poly_1D
Numerically:  4
Analytically: 4

Gauss::Scheme::Quadrangle_5D_7P
Quadrangle_Poly_2D
Numerically:  6.666666666666667
Analytically: 6.666666666666667

Gauss::Scheme::Quadrangle_5D_7P
Quadrangle_Poly_3D
Numerically:  6.666666666666667
Analytically: 6.666666666666667

Gauss::Scheme::Quadrangle_5D_7P
Quadrangle_Poly_4D
Numerically:  8.7111111111111121
Analytically: 8.7111111111111104

Gauss::Scheme::Quadrangle_5D_7P
Quadrangle_Poly_5D
Numerically:  8.7111111111111121
Analytically: 8.7111111111111104

Gauss::Scheme::Quadrangle_7D_12P
Quadrangle_Poly_1D
Numerically:  4
Analytically: 4

Gauss::Scheme::Quadrangle_7D_12P
Quadrangle_Poly_2D
Numerically:  6.666666666666667
Analytically: 6.666666666666667

Gauss::Scheme::Quadrangle_7D_12P
Quadrangle_Poly_3D
Numerically:  6.6666666666666661
Analytically: 6.666666666666667

Gauss::Scheme::Quadrangle_7D_12P
Quadrangle_Poly_4D
Numerically:  8.7111111111111104
Analytically: 8.7111111111111104

Gauss::Scheme::Quadrangle_7D_12P
Quadrangle_Poly_5D
Numerically:  8.7111111111111121
Analytically: 8.7111111111111104

Gauss::Scheme::Quadrangle_7D_12P
Quadrangle_Poly_6D
Numerically:  10.387301587301589
Analytically: 10.387301587301588

Gauss::Scheme::Quadrangle_7D_12P
Quadrangle_Poly_7D
Numerically:  10.387301587301586
Analytically: 10.387301587301588

Gauss::Scheme::Quadrangle_9D_17P
Quadrangle_Poly_1D
Numerically:  4
Analytically: 4

Gauss::Scheme::Quadrangle_9D_17P
Quadrangle_Poly_2D
Numerically:  6.666666666666667
Analytically: 6.666666666666667

Gauss::Scheme::Quadrangle_9D_17P
Quadrangle_Poly_3D
Numerically:  6.666666666666667
Analytically: 6.666666666666667

Gauss::Scheme::Quadrangle_9D_17P
Quadrangle_Poly_4D
Numerically:  8.7111111111111104
Analytically: 8.7111111111111104

Gauss::Scheme::Quadrangle_9D_17P
Quadrangle_Poly_5D
Numerically:  8.7111111111111121
Analytically: 8.7111111111111104

Gauss::Scheme::Quadrangle_9D_17P
Quadrangle_Poly_6D
Numerically:  10.387301587301588
Analytically: 10.387301587301588

Gauss::Scheme::Quadrangle_9D_17P
Quadrangle_Poly_7D
Numerically:  10.387301587301588
Analytically: 10.387301587301588

Gauss::Scheme::Quadrangle_9D_17P
Quadrangle_Poly_8D
Numerically:  11.817142857142859
Analytically: 11.817142857142857

Gauss::Scheme::Quadrangle_9D_17P
Quadrangle_Poly_9D
Numerically:  11.817142857142857
Analytically: 11.817142857142857

Стиль выбран по результатам теста из Поста #351 - без потери производительности, в т.ч. при единственной точке.

[румата]

Сообщение от Нубий-IV Оказывается, гаусс-то засекречен!

Кем он засекречен? https://ru.wikipedia.org/wiki/%D0%A1...BC%D1%83%D0%BB

----- добавлено через -----

Сообщение от Нубий-IV В книгах по интегрированию либо теория с примерами для низких степеней, либо числовые коэффициенты с 5-15 знаками, что недотягивает до точности double.

По ссылке весовые к-ты даны в натуральных дробях. Точности хоть отбавляй

[Нубий-IV]

Еще неприятность. 32-битная программа не только систему уравнений в памяти не держит, но и конечные элементы с результатами расчета.

Код:

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

 
#include <iostream>

const size_t NUM_NODES = 200000;
const size_t NUM_ELEMENTS = 300000;
const size_t NUM_ELEMENT_NODES = 4;
const size_t NUM_NODE_VARS = 6;
const size_t NUM_LOADS = 20;

struct Node {
	double x;
	double y;
	double z;
};

struct Element {
	double K[NUM_NODE_VARS * NUM_ELEMENT_NODES][NUM_NODE_VARS * NUM_ELEMENT_NODES];
};

struct Variable {
	double x;
	double y;
	double z;
	double Ux;
	double Uy;
	double Uz;
};

struct Result {
	double Nx;
	double Ny;
	double Txy;
	double Mx;
	double My;
	double Mxy;

	double Sxx;
	double Syy;
	double Szz;
	double Sxy;
	double Syz;
	double Szx;

	double Asxt;
	double Asxb;
	double Asyt;
	double Asyb;
};

int main()
{
	Node** nodes = new Node*[NUM_NODES];
	for(size_t i = 0; i < NUM_NODES; i++)
		nodes[i] = new Node;
	std::cout << "nodes = " << nodes << std::endl;
	//std::cin.get();

	Element** elements = new Element*[NUM_ELEMENTS];
	for(size_t i = 0; i < NUM_ELEMENTS; i++)
		elements[i] = new Element;
	std::cout << "elements = " << elements << std::endl;
	//std::cin.get();

	Variable** variables = new Variable*[NUM_NODES * NUM_LOADS];
	for(size_t i = 0; i < NUM_NODES * NUM_LOADS; i++)
		variables[i] = new Variable;
	std::cout << "variables = " << variables << std::endl;
	//std::cin.get();

	Result** results = new Result*[NUM_NODES * NUM_LOADS];
	for(size_t i = 0; i < NUM_NODES * NUM_LOADS; i++)
		results[i] = new Result;
	std::cout << "results = " << results << std::endl;
	//std::cin.get();

	/*** а так даже не компилируется при 500000 элементах
	Node* nodes = new Node[NUM_NODES];
	std::cout << "nodes = " << nodes << std::endl;

	Element* elements = new Element[NUM_ELEMENTS];
	std::cout << "elements = " << elements << std::endl;

	Variable* variables = new Variable[NUM_NODES * NUM_LOADS];
	std::cout << "variables = " << variables << std::endl;

	Result* results = new Result[NUM_NODES * NUM_LOADS];
	std::cout << "results = " << results << std::endl;
	*/

	return 0;
}

200 тыс узлов - слишком много. Либо x32 вообще отпадает, либо нужен свой диспетчер данных с подкачкой.

Сообщение от румата По ссылке весовые к-ты даны в натуральных дробях.

Только для нескольких степеней. Где там хотя бы 4 степень для треугольника? Найти что-то выше третьей - пятой степени оказалось проблемой. А если у нас будет нелинейность, повышенные степени могут понадобиться.

[румата]

Сообщение от Нубий-IV Где там хотя бы 4 степень для треугольника?

А для чего треугольнику порядок выше третьего, когда функция формы трегольника полином степени не выше третьей? Этого вполне достаточно для решения как геометрически так и физически нелинейных задач с высокой точностью. Но если вам этого мало - пожалуйста вычисляйте нужное вам количество корней из соответсвующей степени полиномов Лежандра.

[Нубий-IV]

Сообщение от румата А для чего треугольнику порядок выше третьего

Я не писал своих элементов, и не знаю, какие степени будут нужны. Но в исходниках GetDP для треугольника последняя схема - 16-точечная; кому-то оно надо было.

Сообщение от румата вычисляйте нужное вам количество корней

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

Это все. Как по этому описанию получить координаты точек? Ссылка [5,a] - диссертация, нет в свободном доступе. Таблица 15.2 - точки и веса точностью в 7знаков.

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

Кстати, 4-точечная формула из Википедии - единственная возможная, или есть версии без отрицательных коэффициентов?

[румата]

Сообщение от Нубий-IV Кстати, 4-точечная формула из Википедии - единственная возможная, или есть версии без отрицательных коэффициентов?

Я семиточечную применял, тоже вполне рабочий вариант.

Сообщение от Нубий-IV Вопрос-то в том, почему на второй сотне лет численного интегрирования нет полного справочника с формулами и значениями машинной точности.

Может такие справочники и есть где-то, но кому нужны эти сложности, когда для треугольного КЭ вполне достаточно интегрирования 3-го порядка.