mkse

.gitattributes

@@ -0,0 +1,15 @@
+* text=auto
+
+*.cpp text
+*.h text
+*.txt text
+*.md text
+*.cmake text
+*.sln text eol=crlf
+*.vcxproj text eol=crlf
+*.filters text eol=crlf
+
+*.png binary
+*.jpg binary
+*.jpeg binary
+*.pdf binary

.gitignore

@@ -0,0 +1,34 @@
+# Build trees and compiler output
+/build/
+/out/
+/x64/
+/x86/
+/Debug/
+/Release/
+*.exe
+*.ilk
+*.iobj
+*.ipdb
+*.obj
+*.pdb
+*.tlog
+
+# IDE and user-specific files
+/.vs/
+/.vscode/
+*.suo
+*.user
+*.VC.db
+*.VC.opendb
+
+# Generated numerical results
+/res/
+/results/
+
+# Temporary files
+*.log
+*.tmp
+*.bak
+*~
+__pycache__/
+*.pyc

CMakeLists.txt

@@ -0,0 +1,29 @@
+cmake_minimum_required(VERSION 3.20)
+
+project(
+    mkse_elasticity
+    VERSION 1.0.0
+    DESCRIPTION "Finite superelement solver for elastic contact problems"
+    LANGUAGES CXX
+)
+
+add_executable(
+    mkse-elasticity
+    src/main.cpp
+    src/FEM.cpp
+    src/FSEM.cpp
+    src/LinearAlgebra.cpp
+)
+
+target_include_directories(
+    mkse-elasticity
+    PRIVATE "${CMAKE_CURRENT_SOURCE_DIR}/include"
+)
+
+target_compile_features(mkse-elasticity PRIVATE cxx_std_20)
+
+if(MSVC)
+    target_compile_options(mkse-elasticity PRIVATE /W4 /permissive-)
+else()
+    target_compile_options(mkse-elasticity PRIVATE -Wall -Wextra -Wpedantic)
+endif()

include/ContactSlaveNodes.h

@@ -0,0 +1,91 @@
+#pragma once
+
+// Contact discretization based on the passive body's trace nodes.
+
+#include <algorithm>
+#include <cmath>
+#include <stdexcept>
+
+#include "Elasticity.h"
+
+namespace contact_slave_nodes {
+
+constexpr double kContactEps = 1e-12;
+
+inline bool almost_equal(double lhs, double rhs) {
+	return std::fabs(lhs - rhs) < kContactEps;
+}
+
+inline std::vector<double> collect_overlap_nodes(
+	const std::vector<double>& bottom_x,
+	const std::vector<double>& top_x,
+	double contact_left,
+	double contact_right) {
+
+	std::vector<double> nodes;
+	nodes.reserve(bottom_x.size() + top_x.size() + 2);
+	nodes.push_back(contact_left);
+	nodes.push_back(contact_right);
+
+	auto append_inside = [&](const std::vector<double>& x_nodes) {
+		for (double x : x_nodes)
+			if (x >= contact_left - kContactEps && x <= contact_right + kContactEps)
+				nodes.push_back(x);
+	};
+
+	append_inside(bottom_x);
+	append_inside(top_x);
+
+	std::sort(nodes.begin(), nodes.end());
+	nodes.erase(std::unique(nodes.begin(), nodes.end(), almost_equal), nodes.end());
+	return nodes;
+}
+
+inline std::vector<double> collect_active_lambda_nodes(
+	const std::vector<double>& passive_x,
+	double contact_left,
+	double contact_right) {
+
+	std::vector<double> active_nodes;
+	active_nodes.reserve(passive_x.size());
+
+	for (size_t i = 0; i < passive_x.size(); ++i) {
+		const double support_left = (i == 0) ? passive_x[i] : passive_x[i - 1];
+		const double support_right = (i + 1 == passive_x.size()) ? passive_x[i] : passive_x[i + 1];
+
+		if (support_right > contact_left + kContactEps &&
+			support_left < contact_right - kContactEps) {
+			active_nodes.push_back(passive_x[i]);
+		}
+	}
+
+	if (active_nodes.size() < 2)
+		throw std::runtime_error("Not enough passive contact nodes for Lagrange multipliers.");
+
+	return active_nodes;
+}
+
+inline ContactDiscretization build(
+	const std::vector<double>& bottom_x,
+	const std::vector<double>& top_x,
+	double contact_left,
+	double contact_right,
+	ContactSlaveBody passive_body) {
+
+	const std::vector<double>& passive_x =
+		(passive_body == ContactSlaveBody::Bottom) ? bottom_x : top_x;
+
+	ContactDiscretization discretization;
+	discretization.lambda_nodes = collect_active_lambda_nodes(
+		passive_x, contact_left, contact_right);
+
+	const std::vector<double> integration_nodes = collect_overlap_nodes(
+		bottom_x, top_x, contact_left, contact_right);
+	discretization.mortar_elements.reserve(integration_nodes.size() - 1);
+	for (size_t seg = 0; seg + 1 < integration_nodes.size(); ++seg)
+		discretization.mortar_elements.push_back({ integration_nodes[seg], integration_nodes[seg + 1] });
+
+	return discretization;
+}
+
+} // namespace contact_slave_nodes

include/ContactUniformLambdaPartition.h

@@ -0,0 +1,81 @@
+#pragma once
+
+// Contact discretization on an independent uniform multiplier grid.
+
+#include <algorithm>
+#include <cmath>
+
+#include "Elasticity.h"
+
+namespace contact_uniform_lambda_partition {
+
+constexpr double kContactEps = 1e-12;
+
+inline size_t count_contact_nodes(
+	const std::vector<double>& x_nodes,
+	double contact_left,
+	double contact_right) {
+
+	size_t count = 0;
+	for (double x : x_nodes)
+		if (x >= contact_left - kContactEps && x <= contact_right + kContactEps)
+			++count;
+
+	return count;
+}
+
+inline size_t default_lambda_node_count(
+	const std::vector<double>& bottom_x,
+	const std::vector<double>& top_x,
+	double contact_left,
+	double contact_right) {
+
+	return std::max<size_t>(
+		2,
+		std::max(
+			count_contact_nodes(bottom_x, contact_left, contact_right),
+			count_contact_nodes(top_x, contact_left, contact_right)));
+}
+
+inline std::vector<double> build_uniform_lambda_nodes(
+	double contact_left,
+	double contact_right,
+	size_t lambda_node_count) {
+
+	std::vector<double> lambda_nodes(lambda_node_count);
+	const double step = (contact_right - contact_left) / (lambda_node_count - 1);
+
+	for (size_t i = 0; i < lambda_node_count; ++i)
+		lambda_nodes[i] = contact_left + i * step;
+
+	lambda_nodes.front() = contact_left;
+	lambda_nodes.back() = contact_right;
+	return lambda_nodes;
+}
+
+inline ContactDiscretization build(
+	const std::vector<double>& bottom_x,
+	const std::vector<double>& top_x,
+	double contact_left,
+	double contact_right,
+	size_t lambda_node_count) {
+
+	if (lambda_node_count == 0)
+		lambda_node_count = default_lambda_node_count(
+			bottom_x, top_x, contact_left, contact_right);
+
+	ContactDiscretization discretization;
+	discretization.lambda_nodes = build_uniform_lambda_nodes(
+		contact_left, contact_right, std::max<size_t>(2, lambda_node_count));
+
+	discretization.mortar_elements.reserve(discretization.lambda_nodes.size() - 1);
+	for (size_t seg = 0; seg + 1 < discretization.lambda_nodes.size(); ++seg)
+		discretization.mortar_elements.push_back({
+			discretization.lambda_nodes[seg],
+			discretization.lambda_nodes[seg + 1]
+		});
+
+	return discretization;
+}
+
+} // namespace contact_uniform_lambda_partition

include/ContactUniformUnionPartition.h

@@ -0,0 +1,73 @@
+#pragma once
+
+// Contact discretization on a uniform refinement of the combined traces.
+
+#include <algorithm>
+#include <cmath>
+
+#include "Elasticity.h"
+#include "ContactUniformLambdaPartition.h"
+
+namespace contact_uniform_union_partition {
+
+constexpr double kContactEps = 1e-12;
+
+inline bool almost_equal(double lhs, double rhs) {
+	return std::fabs(lhs - rhs) < kContactEps;
+}
+
+inline std::vector<MortarElement> build_mortar_elements(
+	const std::vector<double>& bottom_x,
+	const std::vector<double>& top_x,
+	const std::vector<double>& lambda_nodes,
+	double contact_left,
+	double contact_right) {
+
+	std::vector<double> partition;
+	partition.reserve(bottom_x.size() + top_x.size() + lambda_nodes.size() + 2);
+	partition.push_back(contact_left);
+	partition.push_back(contact_right);
+
+	auto append_contact_nodes = [&](const std::vector<double>& nodes) {
+		for (double x : nodes)
+			if (x >= contact_left - kContactEps && x <= contact_right + kContactEps)
+				partition.push_back(x);
+	};
+
+	append_contact_nodes(bottom_x);
+	append_contact_nodes(top_x);
+	append_contact_nodes(lambda_nodes);
+
+	std::sort(partition.begin(), partition.end());
+	partition.erase(std::unique(partition.begin(), partition.end(), almost_equal), partition.end());
+
+	std::vector<MortarElement> mortar_elements;
+	mortar_elements.reserve(partition.size() - 1);
+	for (size_t seg = 0; seg + 1 < partition.size(); ++seg)
+		mortar_elements.push_back({ partition[seg], partition[seg + 1] });
+
+	return mortar_elements;
+}
+
+inline ContactDiscretization build(
+	const std::vector<double>& bottom_x,
+	const std::vector<double>& top_x,
+	double contact_left,
+	double contact_right,
+	size_t lambda_node_count) {
+
+	ContactDiscretization discretization =
+		contact_uniform_lambda_partition::build(
+			bottom_x, top_x, contact_left, contact_right, lambda_node_count);
+
+	discretization.mortar_elements = build_mortar_elements(
+		bottom_x,
+		top_x,
+		discretization.lambda_nodes,
+		contact_left,
+		contact_right);
+
+	return discretization;
+}
+
+} // namespace contact_uniform_union_partition

include/Elasticity.h

@@ -0,0 +1,265 @@
+#pragma once
+
+// Core data structures and solvers for the finite superelement model.
+
+#include <algorithm>
+#include <cmath>
+#include <cstddef>
+#include <functional>
+#include <iostream>
+#include <string>
+#include <tuple>
+#include <utility>
+#include <vector>
+
+struct Point {
+	double x, y;
+
+	friend std::ostream& operator<<(std::ostream& output, const Point& p);
+};
+
+Point operator+(const Point& a, const Point& b);
+Point operator-(const Point& a, const Point& b);
+Point operator*(double a, const Point& p);
+Point operator/(const Point& p, double a);
+
+std::vector<double> operator-(const std::vector<double>& a,
+	const std::vector<double>& b);
+
+using function = std::function<double(const Point&)>;
+using vec_function = std::function<Point(const Point&)>;
+using lambda_func = std::function<size_t(size_t)>;
+
+enum class CoordinateSystem {
+	Cartesian,
+	Axisymmetric
+};
+
+enum class ContactMethod {
+	SlaveNodes,
+	UniformLambdaPartition,
+	UniformUnionPartition
+};
+
+enum class ContactSlaveBody {
+	Bottom,
+	Top
+};
+
+struct ContactOptions {
+	CoordinateSystem coordinate_system = CoordinateSystem::Axisymmetric;
+	ContactMethod method = ContactMethod::UniformUnionPartition;
+	size_t lambda_node_count = 0;
+	ContactSlaveBody slave_body = ContactSlaveBody::Bottom;
+};
+
+std::string to_string(CoordinateSystem coordinate_system);
+std::string to_string(ContactMethod contact_method);
+std::string to_string(ContactSlaveBody slave_body);
+
+struct Triangle {
+	size_t a, b, c;
+};
+
+struct Rectangle {
+	size_t a, b, c, d;
+};
+
+struct MortarElement {
+	double chi_left;
+	double chi_right;
+};
+
+struct ContactDiscretization {
+	std::vector<double> lambda_nodes;
+	std::vector<MortarElement> mortar_elements;
+};
+
+class Matrix {
+private:
+	std::vector<std::vector<double>> matrix;
+
+public:
+	Matrix(size_t n, double a = 0);
+	Matrix(size_t n, size_t m, double a = 0);
+	Matrix(std::vector<std::vector<double>> m) : matrix(m) {};
+	Matrix(Matrix* M) { matrix = M->matrix; };
+
+	std::vector<double>& operator[](size_t i) { return matrix[i]; };
+	const std::vector<double>& operator[](size_t i) const { return matrix[i]; };
+	Matrix operator*(double a);
+	Matrix& operator*=(double a);
+
+	Matrix T() const;
+	size_t size(short axis = 0) const {
+		return axis ? matrix[0].size() : matrix.size();
+	};
+	Matrix dot(const Matrix& m) const;
+	std::vector<double> dot(const std::vector<double>& v) const;
+	void print() const;
+
+	static Matrix eye(size_t n, double a = 1);
+};
+
+std::vector<double> solveGaussFullPivot(
+	const Matrix& A,
+	const std::vector<double>& b,
+	double eps = 1e-12
+);
+
+class FEM {
+	size_t mx, ny;
+	Point left_down;
+	Point right_up;
+	std::vector<Point> points;
+	std::vector<Triangle> triangles;
+	double triangle_area;
+	std::vector<Point> f;
+	std::vector<Point> u;
+	CoordinateSystem coordinate_system;
+
+	std::vector<double> F;
+	Matrix A;
+
+public:
+	FEM(const Point& a, const Point& b, size_t n, size_t m,
+		CoordinateSystem coordinate_system = CoordinateSystem::Axisymmetric);
+
+	Point& get_point(size_t n);
+	const Point& get_point(size_t n) const;
+	Triangle& get_triangle(size_t n);
+	const Triangle& get_triangle(size_t n) const;
+
+	size_t psize() const { return points.size(); };
+	size_t tsize() const { return triangles.size(); };
+	size_t xsize() const { return mx; };
+	size_t ysize() const { return ny; };
+	CoordinateSystem get_coordinate_system() const { return coordinate_system; };
+
+	void print_points() const;
+	void print_triangles() const;
+
+	Point& operator[](size_t n) { return get_point(n); };
+	const Point& operator[](size_t n) const { return get_point(n); };
+
+	Triangle& operator()(size_t n) { return get_triangle(n); };
+	const Triangle& operator()(size_t n) const { return get_triangle(n); };
+
+	void set_boundaries(char side, const vec_function& g);
+
+	void construct_AF(double E, double nu, vec_function body_force);
+	void apply_boundaries();
+	std::vector<Point> solve();
+
+	void clear_AFu();
+
+	void bc2_side(lambda_func j, size_t start, size_t finish, double len,
+		int side, const std::vector<vec_function>& g,
+		std::vector<double>& p_vec);
+
+	void calculate_bc2(const std::vector<size_t>& pos,
+		const std::vector<vec_function>& g, std::vector<double>& p_vec);
+
+	std::pair<Matrix, std::vector<double>> get_AF();
+	void set_AF(const Matrix& A_new, const std::vector<double>& F_new);
+};
+
+Matrix operator*(double a, Matrix m);
+
+std::tuple<Matrix, Matrix> LU_decomposition(const Matrix& m);
+
+std::vector<double> solveLU(const Matrix& L, const Matrix& U,
+	const std::vector<double>& b);
+
+Point zero(const Point& p);
+
+class FSEM {
+	double E;
+	double nu;
+	Point a;
+	Point b;
+	size_t n_side_x;
+	size_t n_side_y;
+	std::vector<Point> nodes;
+	size_t coef_x;
+	size_t coef_y;
+
+	std::vector<std::vector<Point>> basis;
+	std::vector<Point> basis_coefficients;
+
+	Matrix K;
+	std::vector<double> f;
+	CoordinateSystem coordinate_system;
+
+public:
+	FEM fem;
+
+	FSEM(double E, double nu, const Point& a, const Point& b,
+		size_t n_x, size_t n_y, int coef_val_x = 1, int coef_val_y = 1,
+		CoordinateSystem coordinate_system = CoordinateSystem::Axisymmetric);
+
+	Point& get_node(size_t n);
+	const Point& get_node(size_t n) const;
+
+	size_t nsize() const { return nodes.size(); };
+	CoordinateSystem get_coordinate_system() const { return coordinate_system; };
+
+	Point& operator[](size_t n) { return get_node(n); };
+	const Point& operator[](size_t n) const { return get_node(n); };
+
+	void print_nodes() const;
+
+	void construct_basis();
+
+	Matrix matrix_form_basis();
+
+	const Matrix& get_K() const { return K; }
+	const std::vector<double>& get_f() const { return f; }
+	const std::vector<std::vector<Point>>& get_basis() const { return basis; }
+
+	void construct_f_bc2(const std::vector<size_t>& pos,
+		const std::vector<vec_function>& g);
+
+	std::vector<Point> find_answer();
+	std::vector<Point> find_answer(const std::vector<double>& coefs, size_t start = 0);
+
+	void set_bc1(char side, const vec_function& g);
+
+	void calculate_coef_Matrix_bc2(const int finish, const int i,
+		bool cur_pos, bool prev_pos, int& add_B, int& add_C, const int add_basis,
+		Matrix& B, Matrix& C);
+
+	void save_bc1(std::vector<double>& coefs_Dirichle, int& dir_id, const int i);
+
+	void set_bc2(const std::vector<size_t>& pos,
+		const std::vector<vec_function>& g);
+
+	std::vector<size_t> get_side_nodes(char side) const;
+	std::vector<size_t> get_side_fem_nodes(char side) const;
+
+	Point coefficient(int i, Point coef_val) {
+		return std::isnan(basis_coefficients[i].x) ? coef_val : basis_coefficients[i];
+	}
+
+	std::vector<std::pair<size_t, double>> get_known_dofs() const;
+};
+
+double mortar_shape_func(size_t i, const std::vector<double>& s, double cur);
+
+std::vector<double> solve_mortar_contact(
+	FSEM& bottom_body,
+	FSEM& top_body,
+	const std::vector<double>& rhs_bottom,
+	const std::vector<double>& rhs_top,
+	const ContactOptions& options);
+
+std::vector<double> solve_mortar_contact(
+	FSEM& bottom_body,
+	FSEM& top_body,
+	const std::vector<double>& rhs_bottom,
+	const std::vector<double>& rhs_top,
+	size_t lambda_node_count = 0);
+
+std::vector<double> solveWithLU(const Matrix& A,
+	const std::vector<double>& b,
+	double eps = 1e-15);

include/FEMAxisymmetric.h

@@ -0,0 +1,115 @@
+#pragma once
+
+// Axisymmetric finite-element integration kernels.
+
+#include <array>
+#include <cmath>
+
+#include "Elasticity.h"
+
+namespace fem_axisymmetric {
+
+constexpr double kTwoPi = 6.28318530717958647692;
+
+struct TriangleQuadraturePoint {
+	double weight;
+	std::array<double, 3> phi;
+};
+
+inline const std::array<TriangleQuadraturePoint, 3> kTriangleQuadrature = { {
+	{ 1.0 / 3.0, { 1.0 / 6.0, 1.0 / 6.0, 2.0 / 3.0 } },
+	{ 1.0 / 3.0, { 1.0 / 6.0, 2.0 / 3.0, 1.0 / 6.0 } },
+	{ 1.0 / 3.0, { 2.0 / 3.0, 1.0 / 6.0, 1.0 / 6.0 } }
+} };
+
+inline double signed_double_area(const Point& p1, const Point& p2, const Point& p3) {
+	return (p2.x - p1.x) * (p3.y - p1.y) - (p3.x - p1.x) * (p2.y - p1.y);
+}
+
+inline void assemble(
+	Matrix& A,
+	std::vector<double>& F,
+	const std::vector<Point>& points,
+	const std::vector<Triangle>& triangles,
+	double,
+	double E,
+	double nu,
+	const vec_function& f) {
+
+	const double lambda = (E * nu) / ((1.0 + nu) * (1.0 - 2.0 * nu));
+	const double mu = E / (2.0 * (1.0 + nu));
+	const Matrix C({
+		{ lambda + 2.0 * mu, lambda, lambda, 0.0 },
+		{ lambda, lambda + 2.0 * mu, lambda, 0.0 },
+		{ lambda, lambda, lambda + 2.0 * mu, 0.0 },
+		{ 0.0, 0.0, 0.0, mu }
+	});
+
+	for (const Triangle& T : triangles) {
+		const Point& p1 = points[T.a];
+		const Point& p2 = points[T.b];
+		const Point& p3 = points[T.c];
+		const double two_area = signed_double_area(p1, p2, p3);
+		const double area = 0.5 * std::fabs(two_area);
+
+		const Matrix grad_phi({
+			{ (p2.y - p3.y) / two_area, (p3.x - p2.x) / two_area },
+			{ (p3.y - p1.y) / two_area, (p1.x - p3.x) / two_area },
+			{ (p1.y - p2.y) / two_area, (p2.x - p1.x) / two_area }
+		});
+
+		Matrix Ae(6ull);
+		std::vector<double> Fe(6, 0.0);
+		for (const auto& qp : kTriangleQuadrature) {
+			const Point q = qp.phi[0] * p1 + qp.phi[1] * p2 + qp.phi[2] * p3;
+			const double r_q = q.x;
+
+			Matrix B(4ull, 6ull);
+			for (size_t p = 0; p < 3; ++p) {
+				const double dphi_dr = grad_phi[p][0];
+				const double dphi_dz = grad_phi[p][1];
+				const double phi = qp.phi[p];
+
+				B[0][2 * p] = dphi_dr;
+				B[1][2 * p + 1] = dphi_dz;
+				B[2][2 * p] = phi / r_q;
+				B[3][2 * p] = dphi_dz;
+				B[3][2 * p + 1] = dphi_dr;
+			}
+
+			const Matrix stiffness_q = B.T().dot(C).dot(B);
+			const double weight = kTwoPi * area * qp.weight * r_q;
+			for (size_t i = 0; i < 6; ++i)
+				for (size_t j = 0; j < 6; ++j)
+					Ae[i][j] += weight * stiffness_q[i][j];
+
+			const Point body_force = f(q);
+			for (size_t p = 0; p < 3; ++p) {
+				const double shape_value = qp.phi[p];
+				Fe[2 * p] += weight * shape_value * body_force.x;
+				Fe[2 * p + 1] += weight * shape_value * body_force.y;
+			}
+		}
+
+		auto dof = [T](size_t i) -> size_t {
+			if (i == 0) return 2 * T.a;
+			if (i == 1) return 2 * T.a + 1;
+			if (i == 2) return 2 * T.b;
+			if (i == 3) return 2 * T.b + 1;
+			if (i == 4) return 2 * T.c;
+			return 2 * T.c + 1;
+		};
+
+		for (size_t i = 0; i < 6; ++i) {
+			F[dof(i)] += Fe[i];
+			for (size_t j = 0; j < 6; ++j)
+				A[dof(i)][dof(j)] += Ae[i][j];
+		}
+	}
+}
+
+inline double boundary_segment_weight(const Point& midpoint, double len) {
+	return 0.5 * kTwoPi * midpoint.x * len;
+}
+
+} // namespace fem_axisymmetric

include/FEMCartesian.h

@@ -0,0 +1,76 @@
+#pragma once
+
+// Cartesian finite-element integration kernels.
+
+#include "Elasticity.h"
+
+namespace fem_cartesian {
+
+inline void assemble(
+	Matrix& A,
+	std::vector<double>& F,
+	const std::vector<Point>& points,
+	const std::vector<Triangle>& triangles,
+	double triangle_area,
+	double E,
+	double nu,
+	const vec_function& f) {
+
+	const double lambda = (E * nu) / ((1.0 + nu) * (1.0 - 2.0 * nu));
+	const double mu = E / (2.0 * (1.0 + nu));
+	const Matrix C({
+		{ lambda + 2.0 * mu, lambda, 0.0 },
+		{ lambda, lambda + 2.0 * mu, 0.0 },
+		{ 0.0, 0.0, mu }
+	});
+
+	for (const Triangle& T : triangles) {
+		const Point& p1 = points[T.a];
+		const Point& p2 = points[T.b];
+		const Point& p3 = points[T.c];
+		const Matrix grad_phi({
+			{ (p2.y - p3.y) / (2.0 * triangle_area), (p3.x - p2.x) / (2.0 * triangle_area) },
+			{ (p3.y - p1.y) / (2.0 * triangle_area), (p1.x - p3.x) / (2.0 * triangle_area) },
+			{ (p1.y - p2.y) / (2.0 * triangle_area), (p2.x - p1.x) / (2.0 * triangle_area) }
+		});
+
+		Matrix R(3, 6ull);
+		for (size_t p = 0; p < 3; ++p) {
+			R[0][2 * p] = grad_phi[p][0];
+			R[1][2 * p + 1] = grad_phi[p][1];
+			R[2][2 * p] = grad_phi[p][1];
+			R[2][2 * p + 1] = grad_phi[p][0];
+		}
+
+		Matrix Ae = triangle_area * R.T().dot(C).dot(R);
+
+		std::vector<double> Fe(6, 0.0);
+		const Point center = (p1 + p2 + p3) / 3.0;
+		const Point body_force = f(center);
+		for (size_t p = 0; p < 3; ++p) {
+			Fe[2 * p] = body_force.x * triangle_area / 3.0;
+			Fe[2 * p + 1] = body_force.y * triangle_area / 3.0;
+		}
+
+		auto dof = [T](size_t i) -> size_t {
+			if (i == 0) return 2 * T.a;
+			if (i == 1) return 2 * T.a + 1;
+			if (i == 2) return 2 * T.b;
+			if (i == 3) return 2 * T.b + 1;
+			if (i == 4) return 2 * T.c;
+			return 2 * T.c + 1;
+		};
+
+		for (size_t i = 0; i < 6; ++i) {
+			F[dof(i)] += Fe[i];
+			for (size_t j = 0; j < 6; ++j)
+				A[dof(i)][dof(j)] += Ae[i][j];
+		}
+	}
+}
+
+inline double boundary_segment_weight(const Point&, double len) {
+	return 0.5 * len;
+}
+
+} // namespace fem_cartesian

include/TestCases.h

@@ -0,0 +1,114 @@
+#pragma once
+
+// Manufactured solutions used to validate the numerical model.
+
+#include <array>
+#include <cmath>
+#include <functional>
+#include <map>
+#include <string>
+
+#include "Elasticity.h"
+
+struct MeshSize {
+	size_t bottom_x;
+	size_t bottom_y;
+	size_t top_x;
+	size_t top_y;
+	size_t lambda_nodes;
+};
+
+struct TestCase {
+	std::function<vec_function(double, double)> exact_solution;
+	double E;
+	double nu;
+	Point bottom_a;
+	Point bottom_b;
+	Point top_a;
+	Point top_b;
+	MeshSize mesh;
+	std::array<char, 3> bottom_dirichlet_sides;
+	std::array<char, 3> top_dirichlet_sides;
+};
+
+const std::map<std::string, TestCase> TESTS = {
+	{
+		"axisymmetric_inverse_r",
+		TestCase{
+			[](double, double) {
+				return [](const Point& p) {
+					return Point{ 5.0 / p.x, 0.0 };
+				};
+			},
+			21e10,
+			0.3,
+			Point{ 1.0, 0.0 },
+			Point{ 3.0, 0.5 },
+			Point{ 1.0, 0.5 },
+			Point{ 3.0, 3.0 },
+			MeshSize{ 18, 18, 18, 18, 18 },
+			std::array<char, 3>{ 'W', 'E', 'S' },
+			std::array<char, 3>{ 'W', 'N', 'E' }
+		}
+	},
+	{
+		"axisymmetric_linear",
+		TestCase{
+			[](double, double) {
+				return [](const Point& p) {
+					return Point{ 2.0 * p.x + 5.0 / p.x, 34.0 * p.y };
+				};
+			},
+			21e10,
+			0.3,
+			Point{ 1.0, 0.0 },
+			Point{ 3.0, 0.5 },
+			Point{ 1.0, 0.5 },
+			Point{ 3.0, 3.0 },
+			MeshSize{ 10, 10, 10, 10, 10 },
+			std::array<char, 3>{ 'W', 'E', 'S' },
+			std::array<char, 3>{ 'W', 'N', 'E' }
+		}
+	},
+	{
+		"cartesian_exp",
+		TestCase{
+			[](double, double) {
+				return [](const Point& p) {
+					return Point{
+						std::exp(p.x) * std::cos(p.y - 0.5),
+						-std::exp(p.x) * std::sin(p.y - 0.5)
+					};
+				};
+			},
+			21e10,
+			0.3,
+			Point{ 0.0, 0.0 },
+			Point{ 1.0, 0.5 },
+			Point{ 0.0, 0.5 },
+			Point{ 1.0, 1.0 },
+			MeshSize{ 6, 6, 10, 10, 10 },
+			std::array<char, 3>{ 'W', 'E', 'S' },
+			std::array<char, 3>{ 'W', 'N', 'E' }
+		}
+	},
+	{
+		"cartesian_linear",
+		TestCase{
+			[](double, double) {
+				return [](const Point& p) {
+					return Point{ -21.0 * p.x, 13.0 * p.y };
+				};
+			},
+			21e10,
+			0.3,
+			Point{ 0.0, 0.0 },
+			Point{ 3.0, 1.0 },
+			Point{ 0.0, 1.0 },
+			Point{ 2.0, 4.0 },
+			MeshSize{ 5, 5, 10, 10, 10 },
+			std::array<char, 3>{ 'W', 'E', 'S' },
+			std::array<char, 3>{ 'W', 'N', 'E' }
+		}
+	}
+};

input.txt

@@ -0,0 +1,29 @@
+# coordinate: a/axisymmetric or c/cartesian
+coordinate = a
+
+# contact: s/slave, d/uniform_lambda, u/uniform_union
+contact = u
+
+# passive body is used only for contact=slave
+passive = top
+
+[axisymmetric_inverse_r]
+bottom_x = 5
+bottom_y = 5
+top_x = 9
+top_y = 9
+lambda = 5
+
+[axisymmetric_inverse_r]
+bottom_x = 5
+bottom_y = 5
+top_x = 9
+top_y = 9
+lambda = 7
+
+[axisymmetric_inverse_r]
+bottom_x = 5
+bottom_y = 5
+top_x = 9
+top_y = 9
+lambda = 9

mkse-elasticity.sln

@@ -0,0 +1,31 @@
+
+Microsoft Visual Studio Solution File, Format Version 12.00
+# Visual Studio Version 17
+VisualStudioVersion = 17.4.33403.182
+MinimumVisualStudioVersion = 10.0.40219.1
+Project("{8BC9CEB8-8B4A-11D0-8D11-00A0C91BC942}") = "mkse-elasticity", "mkse-elasticity.vcxproj", "{9E82B28A-0575-4147-BC8D-B9E212C23900}"
+EndProject
+Global
+	GlobalSection(SolutionConfigurationPlatforms) = preSolution
+		Debug|x64 = Debug|x64
+		Debug|x86 = Debug|x86
+		Release|x64 = Release|x64
+		Release|x86 = Release|x86
+	EndGlobalSection
+	GlobalSection(ProjectConfigurationPlatforms) = postSolution
+		{9E82B28A-0575-4147-BC8D-B9E212C23900}.Debug|x64.ActiveCfg = Debug|x64
+		{9E82B28A-0575-4147-BC8D-B9E212C23900}.Debug|x64.Build.0 = Debug|x64
+		{9E82B28A-0575-4147-BC8D-B9E212C23900}.Debug|x86.ActiveCfg = Debug|Win32
+		{9E82B28A-0575-4147-BC8D-B9E212C23900}.Debug|x86.Build.0 = Debug|Win32
+		{9E82B28A-0575-4147-BC8D-B9E212C23900}.Release|x64.ActiveCfg = Release|x64
+		{9E82B28A-0575-4147-BC8D-B9E212C23900}.Release|x64.Build.0 = Release|x64
+		{9E82B28A-0575-4147-BC8D-B9E212C23900}.Release|x86.ActiveCfg = Release|Win32
+		{9E82B28A-0575-4147-BC8D-B9E212C23900}.Release|x86.Build.0 = Release|Win32
+	EndGlobalSection
+	GlobalSection(SolutionProperties) = preSolution
+		HideSolutionNode = FALSE
+	EndGlobalSection
+	GlobalSection(ExtensibilityGlobals) = postSolution
+		SolutionGuid = {34518B04-F18A-4DED-A748-BDDB0ADA38D7}
+	EndGlobalSection
+EndGlobal

mkse-elasticity.vcxproj

@@ -0,0 +1,165 @@
+<?xml version="1.0" encoding="utf-8"?>
+<Project DefaultTargets="Build" xmlns="http://schemas.microsoft.com/developer/msbuild/2003">
+  <ItemGroup Label="ProjectConfigurations">
+    <ProjectConfiguration Include="Debug|Win32">
+      <Configuration>Debug</Configuration>
+      <Platform>Win32</Platform>
+    </ProjectConfiguration>
+    <ProjectConfiguration Include="Release|Win32">
+      <Configuration>Release</Configuration>
+      <Platform>Win32</Platform>
+    </ProjectConfiguration>
+    <ProjectConfiguration Include="Debug|x64">
+      <Configuration>Debug</Configuration>
+      <Platform>x64</Platform>
+    </ProjectConfiguration>
+    <ProjectConfiguration Include="Release|x64">
+      <Configuration>Release</Configuration>
+      <Platform>x64</Platform>
+    </ProjectConfiguration>
+  </ItemGroup>
+  <PropertyGroup Label="Globals">
+    <VCProjectVersion>16.0</VCProjectVersion>
+    <Keyword>Win32Proj</Keyword>
+    <ProjectGuid>{9e82b28a-0575-4147-bc8d-b9e212c23900}</ProjectGuid>
+    <RootNamespace>mkse_elasticity</RootNamespace>
+    <WindowsTargetPlatformVersion>10.0</WindowsTargetPlatformVersion>
+  </PropertyGroup>
+  <Import Project="$(VCTargetsPath)\Microsoft.Cpp.Default.props" />
+  <PropertyGroup Condition="'$(Configuration)|$(Platform)'=='Debug|Win32'" Label="Configuration">
+    <ConfigurationType>Application</ConfigurationType>
+    <UseDebugLibraries>true</UseDebugLibraries>
+    <PlatformToolset>v143</PlatformToolset>
+    <CharacterSet>Unicode</CharacterSet>
+  </PropertyGroup>
+  <PropertyGroup Condition="'$(Configuration)|$(Platform)'=='Release|Win32'" Label="Configuration">
+    <ConfigurationType>Application</ConfigurationType>
+    <UseDebugLibraries>false</UseDebugLibraries>
+    <PlatformToolset>v143</PlatformToolset>
+    <WholeProgramOptimization>true</WholeProgramOptimization>
+    <CharacterSet>Unicode</CharacterSet>
+  </PropertyGroup>
+  <PropertyGroup Condition="'$(Configuration)|$(Platform)'=='Debug|x64'" Label="Configuration">
+    <ConfigurationType>Application</ConfigurationType>
+    <UseDebugLibraries>true</UseDebugLibraries>
+    <PlatformToolset>v143</PlatformToolset>
+    <CharacterSet>Unicode</CharacterSet>
+  </PropertyGroup>
+  <PropertyGroup Condition="'$(Configuration)|$(Platform)'=='Release|x64'" Label="Configuration">
+    <ConfigurationType>Application</ConfigurationType>
+    <UseDebugLibraries>false</UseDebugLibraries>
+    <PlatformToolset>v143</PlatformToolset>
+    <WholeProgramOptimization>true</WholeProgramOptimization>
+    <CharacterSet>Unicode</CharacterSet>
+  </PropertyGroup>
+  <Import Project="$(VCTargetsPath)\Microsoft.Cpp.props" />
+  <ImportGroup Label="ExtensionSettings">
+  </ImportGroup>
+  <ImportGroup Label="Shared">
+  </ImportGroup>
+  <ImportGroup Label="PropertySheets" Condition="'$(Configuration)|$(Platform)'=='Debug|Win32'">
+    <Import Project="$(UserRootDir)\Microsoft.Cpp.$(Platform).user.props" Condition="exists('$(UserRootDir)\Microsoft.Cpp.$(Platform).user.props')" Label="LocalAppDataPlatform" />
+  </ImportGroup>
+  <ImportGroup Label="PropertySheets" Condition="'$(Configuration)|$(Platform)'=='Release|Win32'">
+    <Import Project="$(UserRootDir)\Microsoft.Cpp.$(Platform).user.props" Condition="exists('$(UserRootDir)\Microsoft.Cpp.$(Platform).user.props')" Label="LocalAppDataPlatform" />
+  </ImportGroup>
+  <ImportGroup Label="PropertySheets" Condition="'$(Configuration)|$(Platform)'=='Debug|x64'">
+    <Import Project="$(UserRootDir)\Microsoft.Cpp.$(Platform).user.props" Condition="exists('$(UserRootDir)\Microsoft.Cpp.$(Platform).user.props')" Label="LocalAppDataPlatform" />
+  </ImportGroup>
+  <ImportGroup Label="PropertySheets" Condition="'$(Configuration)|$(Platform)'=='Release|x64'">
+    <Import Project="$(UserRootDir)\Microsoft.Cpp.$(Platform).user.props" Condition="exists('$(UserRootDir)\Microsoft.Cpp.$(Platform).user.props')" Label="LocalAppDataPlatform" />
+  </ImportGroup>
+  <PropertyGroup Label="UserMacros" />
+  <PropertyGroup>
+    <OutDir>$(ProjectDir)build\msbuild\bin\$(Platform)\$(Configuration)\</OutDir>
+    <IntDir>$(ProjectDir)build\msbuild\obj\$(Platform)\$(Configuration)\</IntDir>
+    <TargetName>mkse-elasticity</TargetName>
+  </PropertyGroup>
+  <ItemDefinitionGroup>
+    <ClCompile>
+      <AdditionalIncludeDirectories>$(ProjectDir)include;%(AdditionalIncludeDirectories)</AdditionalIncludeDirectories>
+    </ClCompile>
+  </ItemDefinitionGroup>
+  <ItemDefinitionGroup Condition="'$(Configuration)|$(Platform)'=='Debug|Win32'">
+    <ClCompile>
+      <WarningLevel>Level4</WarningLevel>
+      <SDLCheck>true</SDLCheck>
+      <PreprocessorDefinitions>WIN32;_DEBUG;_CONSOLE;%(PreprocessorDefinitions)</PreprocessorDefinitions>
+      <ConformanceMode>true</ConformanceMode>
+      <LanguageStandard>stdcpp20</LanguageStandard>
+    </ClCompile>
+    <Link>
+      <SubSystem>Console</SubSystem>
+      <GenerateDebugInformation>true</GenerateDebugInformation>
+    </Link>
+  </ItemDefinitionGroup>
+  <ItemDefinitionGroup Condition="'$(Configuration)|$(Platform)'=='Release|Win32'">
+    <ClCompile>
+      <WarningLevel>Level4</WarningLevel>
+      <FunctionLevelLinking>true</FunctionLevelLinking>
+      <IntrinsicFunctions>true</IntrinsicFunctions>
+      <SDLCheck>true</SDLCheck>
+      <PreprocessorDefinitions>WIN32;NDEBUG;_CONSOLE;%(PreprocessorDefinitions)</PreprocessorDefinitions>
+      <ConformanceMode>true</ConformanceMode>
+      <LanguageStandard>stdcpp20</LanguageStandard>
+    </ClCompile>
+    <Link>
+      <SubSystem>Console</SubSystem>
+      <EnableCOMDATFolding>true</EnableCOMDATFolding>
+      <OptimizeReferences>true</OptimizeReferences>
+      <GenerateDebugInformation>true</GenerateDebugInformation>
+    </Link>
+  </ItemDefinitionGroup>
+  <ItemDefinitionGroup Condition="'$(Configuration)|$(Platform)'=='Debug|x64'">
+    <ClCompile>
+      <WarningLevel>Level4</WarningLevel>
+      <SDLCheck>true</SDLCheck>
+      <PreprocessorDefinitions>_DEBUG;_CONSOLE;%(PreprocessorDefinitions)</PreprocessorDefinitions>
+      <ConformanceMode>true</ConformanceMode>
+      <LanguageStandard>stdcpp20</LanguageStandard>
+    </ClCompile>
+    <Link>
+      <SubSystem>Console</SubSystem>
+      <GenerateDebugInformation>true</GenerateDebugInformation>
+    </Link>
+  </ItemDefinitionGroup>
+  <ItemDefinitionGroup Condition="'$(Configuration)|$(Platform)'=='Release|x64'">
+    <ClCompile>
+      <WarningLevel>Level4</WarningLevel>
+      <FunctionLevelLinking>true</FunctionLevelLinking>
+      <IntrinsicFunctions>true</IntrinsicFunctions>
+      <SDLCheck>true</SDLCheck>
+      <PreprocessorDefinitions>NDEBUG;_CONSOLE;%(PreprocessorDefinitions)</PreprocessorDefinitions>
+      <ConformanceMode>true</ConformanceMode>
+      <LanguageStandard>stdcpp20</LanguageStandard>
+    </ClCompile>
+    <Link>
+      <SubSystem>Console</SubSystem>
+      <EnableCOMDATFolding>true</EnableCOMDATFolding>
+      <OptimizeReferences>true</OptimizeReferences>
+      <GenerateDebugInformation>true</GenerateDebugInformation>
+    </Link>
+  </ItemDefinitionGroup>
+  <ItemGroup>
+    <ClCompile Include="src\FEM.cpp" />
+    <ClCompile Include="src\FSEM.cpp" />
+    <ClCompile Include="src\LinearAlgebra.cpp" />
+    <ClCompile Include="src\main.cpp" />
+  </ItemGroup>
+  <ItemGroup>
+    <ClInclude Include="include\ContactSlaveNodes.h" />
+    <ClInclude Include="include\ContactUniformLambdaPartition.h" />
+    <ClInclude Include="include\ContactUniformUnionPartition.h" />
+    <ClInclude Include="include\Elasticity.h" />
+    <ClInclude Include="include\FEMAxisymmetric.h" />
+    <ClInclude Include="include\FEMCartesian.h" />
+    <ClInclude Include="include\TestCases.h" />
+  </ItemGroup>
+  <ItemGroup>
+    <None Include="input.txt" />
+    <None Include="README.md" />
+  </ItemGroup>
+  <Import Project="$(VCTargetsPath)\Microsoft.Cpp.targets" />
+  <ImportGroup Label="ExtensionTargets">
+  </ImportGroup>
+</Project>

mkse-elasticity.vcxproj.filters

@@ -0,0 +1,54 @@
+<?xml version="1.0" encoding="utf-8"?>
+<Project ToolsVersion="4.0" xmlns="http://schemas.microsoft.com/developer/msbuild/2003">
+  <ItemGroup>
+    <Filter Include="Source Files">
+      <UniqueIdentifier>{4FC737F1-C7A5-4376-A066-2A32D752A2FF}</UniqueIdentifier>
+      <Extensions>cpp;c;cc;cxx;c++;cppm;ixx;def;odl;idl;hpj;bat;asm;asmx</Extensions>
+    </Filter>
+    <Filter Include="Header Files">
+      <UniqueIdentifier>{93995380-89BD-4b04-88EB-625FBE52EBFB}</UniqueIdentifier>
+      <Extensions>h;hh;hpp;hxx;h++;hm;inl;inc;ipp;xsd</Extensions>
+    </Filter>
+    <Filter Include="Resource Files">
+      <UniqueIdentifier>{67DA6AB6-F800-4c08-8B7A-83BB121AAD01}</UniqueIdentifier>
+      <Extensions>rc;ico;cur;bmp;dlg;rc2;rct;bin;rgs;gif;jpg;jpeg;jpe;resx;tiff;tif;png;wav;mfcribbon-ms</Extensions>
+    </Filter>
+  </ItemGroup>
+  <ItemGroup>
+    <ClCompile Include="src\main.cpp">
+      <Filter>Source Files</Filter>
+    </ClCompile>
+    <ClCompile Include="src\FEM.cpp">
+      <Filter>Source Files</Filter>
+    </ClCompile>
+    <ClCompile Include="src\LinearAlgebra.cpp">
+      <Filter>Source Files</Filter>
+    </ClCompile>
+    <ClCompile Include="src\FSEM.cpp">
+      <Filter>Source Files</Filter>
+    </ClCompile>
+  </ItemGroup>
+  <ItemGroup>
+    <ClInclude Include="include\ContactSlaveNodes.h">
+      <Filter>Header Files</Filter>
+    </ClInclude>
+    <ClInclude Include="include\ContactUniformLambdaPartition.h">
+      <Filter>Header Files</Filter>
+    </ClInclude>
+    <ClInclude Include="include\ContactUniformUnionPartition.h">
+      <Filter>Header Files</Filter>
+    </ClInclude>
+    <ClInclude Include="include\FEMAxisymmetric.h">
+      <Filter>Header Files</Filter>
+    </ClInclude>
+    <ClInclude Include="include\FEMCartesian.h">
+      <Filter>Header Files</Filter>
+    </ClInclude>
+    <ClInclude Include="include\Elasticity.h">
+      <Filter>Header Files</Filter>
+    </ClInclude>
+    <ClInclude Include="include\TestCases.h">
+      <Filter>Header Files</Filter>
+    </ClInclude>
+  </ItemGroup>
+</Project>

src/FEM.cpp

@@ -0,0 +1,274 @@
+#include <iostream>
+#include <stdexcept>
+#include <array>
+#include <cmath>
+#include "Elasticity.h"
+#include "FEMAxisymmetric.h"
+#include "FEMCartesian.h"
+
+namespace {
+
+	constexpr double kAxisymmetricTwoPi = 6.28318530717958647692;
+	constexpr double kAxisymmetricRadiusEps = 1e-14;
+
+	struct TriangleQuadraturePoint {
+		double weight;
+		std::array<double, 3> phi;
+	};
+
+	const std::array<TriangleQuadraturePoint, 3> kTriangleQuadrature = { {
+		{ 1.0 / 3.0, { 1.0 / 6.0, 1.0 / 6.0, 2.0 / 3.0 } },
+		{ 1.0 / 3.0, { 1.0 / 6.0, 2.0 / 3.0, 1.0 / 6.0 } },
+		{ 1.0 / 3.0, { 2.0 / 3.0, 1.0 / 6.0, 1.0 / 6.0 } }
+	} };
+
+	double signed_double_area(const Point& p1, const Point& p2, const Point& p3) {
+		return (p2.x - p1.x) * (p3.y - p1.y) - (p3.x - p1.x) * (p2.y - p1.y);
+	}
+
+}
+
+std::ostream& operator<<(std::ostream& output, const Point& p) {
+	output << "{ " << p.x << "; " << p.y << " }";
+	return output;
+}
+
+Point operator+(const Point& a, const Point& b) {
+	return { a.x + b.x, a.y + b.y };
+}
+
+Point operator-(const Point& a, const Point& b) {
+	return { a.x - b.x, a.y - b.y };
+}
+
+Point operator*(const double a, const Point& p) {
+	return { p.x * a, p.y * a };
+}
+
+Point operator/(const Point& p, const double a) {
+	return { p.x / a, p.y / a };
+}
+
+function mul(function a, function b) {
+	return [a, b](const Point& p) { return b(p) * a(p); };
+}
+
+function get_x(vec_function a) {
+	return [a](const Point& p) { return a(p).x; };
+}
+
+function get_y(vec_function a) {
+	return [a](const Point& p) { return a(p).y; };
+}
+
+
+FEM::FEM(const Point& a, const Point& b, size_t n, size_t m,
+	CoordinateSystem coordinate_system)
+	: mx(m), ny(n), A(2 * n * m, 2 * n * m), F(2 * n * m),
+	left_down(a), right_up(b), coordinate_system(coordinate_system) {
+
+	double dx = (b.x - a.x) / (m - 1), dy = (b.y - a.y) / (n - 1);
+	u = std::vector<Point>(m * n, { NAN, NAN });
+	triangle_area = dx * dy / 2;
+	
+	points.resize(m * n);
+	f.resize(m * n);
+	for (size_t i = 0; i != points.size(); ++i) {
+		points[i] = { a.x + (i % m) * dx, a.y + (i / m) * dy };
+	}
+
+	triangles.resize(2 * (m - 1) * (n - 1));
+	for (size_t i = 0; i != triangles.size(); ++i) {
+		size_t x = i % (2 * (m - 1)), y = i / (2 * (m - 1));
+		if (x % 2 == 0) { // "верхний" треугольник
+			triangles[i] = {
+				(y * m) + x / 2,
+				((y + 1) * m) + x / 2,
+				((y + 1) * m) + x / 2 + 1
+			};
+		}
+		else { // "нижний" треугольник
+			triangles[i] = {
+				(y * m) + x / 2,
+				((y + 1) * m) + x / 2 + 1,
+				(y * m) + x / 2 + 1
+			};
+		}
+	}
+}
+
+Point& FEM::get_point(size_t n) {
+	return points[n];
+}
+
+const Point& FEM::get_point(size_t n) const {
+	return points[n];
+}
+
+Triangle& FEM::get_triangle(size_t n) {
+	return triangles[n];
+}
+
+const Triangle& FEM::get_triangle(size_t n) const {
+	return triangles[n];
+}
+
+void FEM::print_points() const {
+	for (size_t i = 0; i < points.size(); ++i) {
+		std::cout << (*this)[i] << "\t";
+	}
+	std::cout << '\n';
+}
+
+void FEM::print_triangles() const {
+	for (size_t i = 0; i < triangles.size(); ++i) {
+		std::cout << "[ " <<
+			(*this)(i).a << ", " <<
+			(*this)(i).b << ", " <<
+			(*this)(i).c << "]\n";
+	}
+}
+
+void FEM::set_boundaries(char side, const vec_function& g) {
+	if (side == 'S') { // Нижняя граница
+		for (size_t j = 0; j != mx; ++j) {
+			Point t = g(points[j]);
+			if(isnan(u[j].x))
+				u[j].x = t.x;
+			if (isnan(u[j].y))
+				u[j].y = t.y;
+		}
+	}
+	else if (side == 'E') { // Правая граница
+		for (size_t i = 0; i != ny; ++i) {
+			Point t = g(points[mx - 1 + i * mx]);
+			if (isnan(u[mx - 1 + i * mx].x))
+				u[mx - 1 + i * mx].x = t.x;
+			if (isnan(u[mx - 1 + i * mx].y))
+				u[mx - 1 + i * mx].y = t.y;
+		}
+	}
+	else if (side == 'N') { // Верхняя граница
+		for (size_t j = 0; j != mx; ++j) {
+			Point t = g(points[mx * (ny - 1) + j]);
+			if (isnan(u[mx * (ny - 1) + j].x))
+				u[mx * (ny - 1) + j].x = t.x;
+			if (isnan(u[mx * (ny - 1) + j].y))
+				u[mx * (ny - 1) + j].y = t.y;
+		}
+	}
+	else { // Левая граница
+		for (size_t i = 0; i != ny; ++i) {
+			Point t = g(points[i * mx]);
+			if (isnan(u[i * mx].x))
+				u[i * mx].x = t.x;
+			if (isnan(u[i * mx].y))
+				u[i * mx].y = t.y;
+		}
+	}
+}
+
+void FEM::construct_AF(double E, double nu, vec_function body_force) {
+	if (coordinate_system == CoordinateSystem::Cartesian) {
+		fem_cartesian::assemble(A, F, points, triangles, triangle_area, E, nu, body_force);
+	}
+	else {
+		fem_axisymmetric::assemble(A, F, points, triangles, triangle_area, E, nu, body_force);
+	}
+}
+
+void FEM::apply_boundaries() {
+	for (size_t i = 0; i < psize(); ++i) {
+		if (!isnan(u[i].x)) {
+			F[2 * i] = u[i].x;
+			for (size_t j = 0; j < A[2 * i].size(); ++j)
+				A[2 * i][j] = 0;
+			A[2 * i][2 * i] = 1;
+		}
+		if (!isnan(u[i].y)) {
+			F[2 * i + 1] = u[i].y;
+			for (size_t j = 0; j < A[2 * i + 1].size(); ++j)
+				A[2 * i + 1][j] = 0;
+			A[2 * i + 1][2 * i + 1] = 1;
+		}
+	}
+}
+
+std::vector<Point> FEM::solve() {
+	auto [L, U] = LU_decomposition(A);
+	std::vector<double> dofs = solveLU(L, U, F);
+
+	std::vector<Point> res(psize());
+	for (size_t i = 0; i < psize(); ++i) {
+		res[i].x = dofs[2 * i];
+		res[i].y = dofs[2 * i + 1];
+	}
+
+	clear_AFu();
+
+	return res;
+}
+
+void FEM::clear_AFu() {
+	for (size_t i = 0; i < 2 * psize(); ++i) {
+		F[i] = 0;
+		if (i % 2 == 0)
+			u[i / 2] = { NAN, NAN };
+		for (size_t j = 0; j < 2 * psize(); ++j)
+			A[i][j] = 0;
+	}
+}
+
+std::pair<Matrix, std::vector<double>> FEM::get_AF() {
+	return { A, F };
+}
+
+void FEM::set_AF(const Matrix& A_new, const std::vector<double>& F_new) {
+	A = A_new;
+	F = F_new;
+}
+
+void FEM::bc2_side(lambda_func j, size_t start, size_t finish, double len,
+	int side, const std::vector<vec_function>& g,
+	std::vector<double>& p_vec) {
+	for (size_t i = start; i < finish; ++i) { // индекс по границе как в мксэ
+		const Point midpoint = (points[j(i)] + points[j(i + 1)]) / 2;
+		const double weight = (coordinate_system == CoordinateSystem::Cartesian)
+			? fem_cartesian::boundary_segment_weight(midpoint, len)
+			: fem_axisymmetric::boundary_segment_weight(midpoint, len);
+		Point integral = weight * g[side](midpoint);
+
+		// текущий узел
+		p_vec[2 * j(i)] += integral.x;
+		p_vec[2 * j(i) + 1] += integral.y;
+
+		p_vec[2 * j(i + 1)] += integral.x;
+		p_vec[2 * j(i + 1) + 1] += integral.y;
+	}
+}
+
+void FEM::calculate_bc2(const std::vector<size_t>& pos,
+	const std::vector<vec_function>& g, std::vector<double>& p_vec) {
+	
+	double len_vert = (right_up.y - left_down.y) / (ny - 1),
+		len_hor = (right_up.x - left_down.x) / (mx - 1);
+
+	if (pos[0])  // слева ГУ 2 рода
+		bc2_side([&](size_t i) { return i * mx; }, 0, ny - 1,
+			len_vert, 0, g, p_vec);
+
+	if (pos[1])  // сверху ГУ 2 рода
+		bc2_side([&](size_t i) { return (mx - 1) * (ny - 1) + i; },
+			ny - 1, ny + mx - 2, len_hor, 1, g, p_vec);
+
+	if (pos[2])  // справа ГУ 2 рода
+		bc2_side([&](size_t i) { return mx * (2 * ny + mx - i - 2) - 1; },
+			ny + mx - 2, 2 * ny + mx - 3, len_vert,
+			2, g, p_vec);
+	
+	if (pos[3])  // снизу ГУ 2 рода
+		bc2_side([&](size_t i) { return 1 + 2 * (ny + mx) - 5 - i; },
+			2 * ny + mx - 3, 2 * (ny + mx) - 4, len_hor,
+			3, g, p_vec);
+
+}

src/LinearAlgebra.cpp

@@ -0,0 +1,328 @@
+#include "Elasticity.h"
+
+#include <stdexcept>
+
+Matrix::Matrix(size_t n, double a) {
+	matrix = std::vector<std::vector<double>>(n, std::vector<double>(n, a));
+}
+
+Matrix::Matrix(size_t n, size_t m, double a) {
+	matrix = std::vector<std::vector<double>>(n, std::vector<double>(m, a));
+}
+
+Matrix Matrix::operator*(double a) {
+	Matrix m = *this;
+
+	for (size_t i = 0; i != matrix.size(); ++i) {
+		for (size_t j = 0; j != matrix[i].size(); ++j) {
+			m[i][j] *= a;
+		}
+	}
+
+	return m;
+}
+
+Matrix operator*(double a, Matrix m) {
+	return m * a;
+}
+
+Matrix& Matrix::operator*=(double a) {
+	for (size_t i = 0; i != matrix.size(); ++i) {
+		for (size_t j = 0; j != matrix[i].size(); ++j) {
+			matrix[i][j] *= a;
+		}
+	}
+
+	return *this;
+}
+
+std::vector<double> operator-(const std::vector<double>& a, 
+	const std::vector<double>& b) {
+	if (a.size() != b.size()) 
+		std::cout << "Wrong vector size for subtraction\n";
+
+	std::vector<double> res(a.size());
+
+	for (int i = 0; i < a.size(); ++i)
+		res[i] = a[i] - b[i];
+
+	return res;
+}
+
+Matrix Matrix::T() const {
+	Matrix m(size(1), size());
+
+	for (size_t i = 0; i != size(1); ++i) {
+		for (size_t j = 0; j != size(); ++j) {
+			m[i][j] = matrix[j][i];
+		}
+	}
+
+	return m;
+}
+
+Matrix Matrix::dot(const Matrix& m) const {
+	Matrix res(size(), m.size(1));
+
+	if (m.size() != size(1))
+		throw std::runtime_error("Wrong matrix size while multiplication!");
+
+	for (size_t i = 0; i != res.size(); ++i) {
+		for (size_t j = 0; j != res.size(1); ++j) {
+
+			for (size_t k = 0; k != m.size(); ++k) {
+				res[i][j] += matrix[i][k] * m[k][j];
+			}
+		}
+	}
+
+	return res;
+}
+
+std::vector<double> Matrix::dot(const std::vector<double>&v) const {
+	if (v.size() != size(1))
+		throw std::runtime_error("Wrong matrix size while multiplication!");
+
+	std::vector<double> res(size());
+
+	for (int i = 0; i < size(); ++i)
+		for (int j = 0; j < size(1); ++j)
+			res[i] += matrix[i][j] * v[j];
+
+	return res;
+}
+
+void Matrix::print() const {
+	for (size_t i = 0; i != size(); ++i) {
+		std::cout << "{";
+		for (size_t j = 0; j != size(1); ++j) {
+			std::cout << matrix[i][j] << ", ";
+		}
+		std::cout << "},\n";
+	}
+}
+
+Matrix Matrix::eye(size_t n, double a) {
+	Matrix m(n, n);
+	for (size_t i = 0; i != m.size(); ++i) {
+		m[i][i] = a;
+	}
+
+	return m;
+}
+
+std::tuple<Matrix, Matrix> LU_decomposition(const Matrix& m) {
+	Matrix U(m.size(), m.size(1));
+	Matrix L = Matrix::eye(m.size());
+
+	for (size_t i = 0; i != m.size(); ++i) {
+		for (size_t j = 0; j != m.size(); ++j) {
+			if (i <= j) {
+				U[i][j] = m[i][j];
+				for (size_t k = 0; k != i; ++k) {
+					U[i][j] -= L[i][k] * U[k][j];
+				}
+			}
+
+			else {
+				L[i][j] = m[i][j];
+				for (size_t k = 0; k != j; ++k) {
+					L[i][j] -= L[i][k] * U[k][j];
+				}
+				L[i][j] /= U[j][j];
+			}
+		}
+	}
+
+	return { L, U };
+}
+
+std::vector<double> solveLU(const Matrix& L, const Matrix& U,
+	const std::vector<double>& b) {
+	std::vector<double> y(L.size());
+	for (size_t i = 0; i != L.size(); ++i) {
+		y[i] = b[i];
+		for (size_t j = 0; j != i; ++j) {
+			y[i] -= L[i][j] * y[j];
+		}
+	}
+
+	std::vector<double> x(L.size());
+	for (size_t i = L.size() - 1; i + 1 != 0; --i) {
+		x[i] = y[i];
+		for (size_t j = i + 1; j != L.size(); ++j) {
+			x[i] -= U[i][j] * x[j];
+		}
+		x[i] /= U[i][i];
+	}
+
+	return x;
+}
+
+std::vector<double> solveGaussFullPivot(const Matrix& A,
+	const std::vector<double>& b,
+	double eps) {
+	
+	const size_t n = A.size();
+	Matrix M(A);
+	std::vector<double> rhs = b;
+
+	std::vector<size_t> col_perm(n);
+
+	for (size_t i = 0; i != n; ++i)
+		col_perm[i] = i;
+
+	for (size_t k = 0; k != n; ++k) {
+		size_t pivot_row = k;
+		size_t pivot_col = k;
+		double pivot_abs = 0;
+
+		for (size_t i = k; i != n; ++i) {
+			for (size_t j = k; j != n; ++j) {
+				double cur = abs(M[i][j]);
+				if (cur > pivot_abs) {
+					pivot_abs = cur;
+					pivot_row = i;
+					pivot_col = j;
+				}
+			}
+		}
+
+		if (pivot_abs < eps)
+			throw std::runtime_error("Cannot solve a singular linear system.");
+
+		if (pivot_row != k) {
+			std::swap(M[pivot_row], M[k]);
+			std::swap(rhs[pivot_row], rhs[k]);
+		}
+
+		if (pivot_col != k) {
+			for (size_t i = 0; i != n; ++i)
+				std::swap(M[i][pivot_col], M[i][k]);
+			std::swap(col_perm[pivot_col], col_perm[k]);
+		}
+
+		for (size_t i = k + 1; i != n; ++i) {
+			double factor = M[i][k] / M[k][k];
+			M[i][k] = 0;
+			for (size_t j = k + 1; j != n; ++j)
+				M[i][j] -= factor * M[k][j];
+			rhs[i] -= factor * rhs[k];
+		}
+	}
+
+	std::vector<double> y(n);
+	for (size_t i = n; i-- > 0;) {
+		double sum = rhs[i];
+		for (size_t j = i + 1; j != n; ++j)
+			sum -= M[i][j] * y[j];
+
+		y[i] = sum / M[i][i];
+	}
+
+	std::vector<double> x(n);
+	for (size_t i = 0; i != n; ++i)
+		x[col_perm[i]] = y[i];
+
+	return x;
+}
+
+bool lu_decomposition_partial_pivot(const Matrix& A,
+	Matrix& L,
+	Matrix& U,
+	std::vector<size_t>& row_perm,
+	double eps = 1e-15) {
+	const size_t n = A.size();
+
+	U = A;
+	L = Matrix::eye(n);
+	row_perm.resize(n);
+	for (size_t i = 0; i < n; ++i) row_perm[i] = i;
+
+	double max_abs = 0.0;
+	for (size_t i = 0; i < n; ++i)
+		for (size_t j = 0; j < n; ++j)
+			max_abs = std::max(max_abs, std::abs(A[i][j]));
+
+	for (size_t k = 0; k < n; ++k) {
+		size_t pivot_row = k;
+		double pivot_abs = 0.0;
+		for (size_t i = k; i < n; ++i) {
+			double cur = std::abs(U[i][k]);
+			if (cur > pivot_abs) { pivot_abs = cur; pivot_row = i; }
+		}
+
+		if (pivot_abs < eps) {
+			double tiny = 1e-12;
+			double tau = tiny * (1.0 + max_abs);
+			for (size_t i = k; i < n; ++i) U[i][i] += tau;
+			pivot_abs = 0.0;
+			pivot_row = k;
+			for (size_t i = k; i < n; ++i) {
+				double cur = std::abs(U[i][k]);
+				if (cur > pivot_abs) { pivot_abs = cur; pivot_row = i; }
+			}
+			if (pivot_abs < eps) return false;
+		}
+
+		if (pivot_row != k) {
+			std::swap(U[pivot_row], U[k]);
+			std::swap(row_perm[pivot_row], row_perm[k]);
+			for (size_t j = 0; j < k; ++j)
+				std::swap(L[pivot_row][j], L[k][j]);
+		}
+
+		double Akk = U[k][k];
+		if (std::abs(Akk) < eps) return false;
+
+		for (size_t i = k + 1; i < n; ++i) {
+			double mult = U[i][k] / Akk;
+			L[i][k] = mult;
+			U[i][k] = 0.0;
+			for (size_t j = k + 1; j < n; ++j)
+				U[i][j] -= mult * U[k][j];
+		}
+	}
+
+	return true;
+}
+
+std::vector<double> solveLU_with_perm(const Matrix& L,
+	const Matrix& U,
+	const std::vector<size_t>& row_perm,
+	const std::vector<double>& b) {
+	const size_t n = L.size();
+	std::vector<double> rhs(n);
+	for (size_t i = 0; i < n; ++i) rhs[i] = b[row_perm[i]];
+
+	std::vector<double> y(n);
+	for (size_t i = 0; i < n; ++i) {
+		double s = rhs[i];
+		for (size_t j = 0; j < i; ++j) s -= L[i][j] * y[j];
+		y[i] = s; 
+	}
+
+	std::vector<double> x(n);
+	for (size_t ii = 0; ii < n; ++ii) {
+		size_t i = n - 1 - ii;
+		double s = y[i];
+		for (size_t j = i + 1; j < n; ++j) s -= U[i][j] * x[j];
+		x[i] = s / U[i][i];
+	}
+
+	return x;
+}
+
+std::vector<double> solveWithLU(const Matrix& A,
+	const std::vector<double>& b,
+	double eps) {
+	const size_t n = A.size();
+
+	Matrix L(n, n), U(n, n);
+	std::vector<size_t> row_perm;
+	bool ok = lu_decomposition_partial_pivot(A, L, U, row_perm, eps);
+	if (!ok) {
+		throw std::runtime_error("LU: pivot ~ 0");
+	}
+	return solveLU_with_perm(L, U, row_perm, b);
+}

src/main.cpp

@@ -0,0 +1,775 @@
+#include "Elasticity.h"
+#include "TestCases.h"
+#include "ContactSlaveNodes.h"
+#include "ContactUniformLambdaPartition.h"
+#include "ContactUniformUnionPartition.h"
+
+#include <algorithm>
+#include <cctype>
+#include <filesystem>
+#include <fstream>
+#include <iomanip>
+#include <iostream>
+#include <optional>
+#include <sstream>
+#include <stdexcept>
+
+namespace {
+
+constexpr double kTraceEps = 1e-12;
+
+struct TestRunConfig {
+	std::string name;
+	std::optional<size_t> bottom_x;
+	std::optional<size_t> bottom_y;
+	std::optional<size_t> top_x;
+	std::optional<size_t> top_y;
+	std::optional<size_t> lambda_nodes;
+};
+
+struct InputConfig {
+	CoordinateSystem coordinate_system = CoordinateSystem::Axisymmetric;
+	ContactMethod contact_method = ContactMethod::UniformUnionPartition;
+	ContactSlaveBody slave_body = ContactSlaveBody::Bottom;
+	std::vector<TestRunConfig> tests;
+};
+
+std::string trim(const std::string& text) {
+	std::string value = text;
+	if (value.size() >= 3 &&
+		static_cast<unsigned char>(value[0]) == 0xEF &&
+		static_cast<unsigned char>(value[1]) == 0xBB &&
+		static_cast<unsigned char>(value[2]) == 0xBF) {
+		value.erase(0, 3);
+	}
+
+	const auto first = value.find_first_not_of(" \t\r\n");
+	if (first == std::string::npos)
+		return {};
+	const auto last = value.find_last_not_of(" \t\r\n");
+	return value.substr(first, last - first + 1);
+}
+
+std::string normalize(std::string value) {
+	value = trim(value);
+	if (value.size() >= 2 &&
+		((value.front() == '"' && value.back() == '"') ||
+			(value.front() == '\'' && value.back() == '\''))) {
+		value = value.substr(1, value.size() - 2);
+	}
+	std::replace(value.begin(), value.end(), '-', '_');
+	for (char& ch : value)
+		ch = static_cast<char>(std::tolower(static_cast<unsigned char>(ch)));
+	return value;
+}
+
+std::vector<std::string> split_list(std::string value) {
+	for (char& ch : value)
+		if (ch == ',' || ch == ';')
+			ch = ' ';
+
+	std::vector<std::string> result;
+	std::istringstream in(value);
+	std::string token;
+	while (in >> token)
+		result.push_back(token);
+	return result;
+}
+
+bool starts_with(const std::string& value, const std::string& prefix) {
+	return value.size() >= prefix.size() &&
+		std::equal(prefix.begin(), prefix.end(), value.begin());
+}
+
+size_t parse_size(const std::string& value, const std::string& key) {
+	size_t parsed = 0;
+	size_t used = 0;
+	parsed = std::stoull(value, &used);
+	if (used != value.size())
+		throw std::runtime_error("Invalid integer value for '" + key + "': " + value);
+	return parsed;
+}
+
+CoordinateSystem parse_coordinate_system(const std::string& value) {
+	const std::string v = normalize(value);
+	if (v == "a" || v == "axi" || v == "axisymmetric") {
+		return CoordinateSystem::Axisymmetric;
+	}
+	if (v == "c" || v == "cartesian" || v == "cart") {
+		return CoordinateSystem::Cartesian;
+	}
+	throw std::runtime_error("Unknown coordinate system: " + value);
+}
+
+ContactMethod parse_contact_method(const std::string& value) {
+	const std::string v = normalize(value);
+	if (v == "s" || v == "slave" || v == "slave_nodes") {
+		return ContactMethod::SlaveNodes;
+	}
+	if (v == "d" || v == "uniform_lambda" ) {
+		return ContactMethod::UniformLambdaPartition;
+	}
+	if (v == "u" || v == "uniform_union") {
+		return ContactMethod::UniformUnionPartition;
+	}
+	throw std::runtime_error("Unknown contact method: " + value);
+}
+
+ContactSlaveBody parse_slave_body(const std::string& value) {
+	const std::string v = normalize(value);
+	if (v == "bottom" || v == "b")
+		return ContactSlaveBody::Bottom;
+	if (v == "top" || v == "t")
+		return ContactSlaveBody::Top;
+	throw std::runtime_error("Unknown passive body: " + value);
+}
+
+TestRunConfig& append_test(InputConfig& config, const std::string& name) {
+	config.tests.push_back(TestRunConfig{ name });
+	return config.tests.back();
+}
+
+void add_test_name(InputConfig& config, const std::string& name) {
+	if (normalize(name) == "all") {
+		for (const auto& [test_name, test] : TESTS) {
+			(void)test;
+			append_test(config, test_name);
+		}
+	}
+	else {
+		append_test(config, name);
+	}
+}
+
+void apply_override(TestRunConfig& run, const std::string& key, const std::string& value) {
+	const std::string k = normalize(key);
+	const std::string v = normalize(value);
+	if (k == "bottom_x")
+		run.bottom_x = parse_size(v, key);
+	else if (k == "bottom_y")
+		run.bottom_y = parse_size(v, key);
+	else if (k == "top_x")
+		run.top_x = parse_size(v, key);
+	else if (k == "top_y")
+		run.top_y = parse_size(v, key);
+	else if (k == "lambda")
+		run.lambda_nodes = parse_size(v, key);
+	else
+		throw std::runtime_error("Unknown per-test parameter: " + key);
+}
+
+void apply_key_value(InputConfig& config, TestRunConfig* current_test,
+	const std::string& key, const std::string& value) {
+
+	const std::string k = normalize(key);
+	if (k == "coordinate") {
+		config.coordinate_system = parse_coordinate_system(value);
+	}
+	else if (k == "contact") {
+		config.contact_method = parse_contact_method(value);
+	}
+	else if (k == "passive") {
+		config.slave_body = parse_slave_body(value);
+	}
+	else if (k == "tests") {
+		for (const std::string& name : split_list(value))
+			add_test_name(config, name);
+	}
+	else if (current_test != nullptr) {
+		apply_override(*current_test, key, value);
+	}
+	else {
+		throw std::runtime_error("Unknown input key outside a test block: " + key);
+	}
+}
+
+void parse_assignment_line(InputConfig& config, TestRunConfig* current_test, const std::string& line) {
+	const size_t eq = line.find('=');
+	const size_t colon = line.find(':');
+	const size_t pos = std::min(
+		eq == std::string::npos ? line.size() : eq,
+		colon == std::string::npos ? line.size() : colon);
+
+	if (pos != line.size()) {
+		apply_key_value(config, current_test, line.substr(0, pos), line.substr(pos + 1));
+		return;
+	}
+
+	std::istringstream in(line);
+	std::string key;
+	std::string value;
+	in >> key >> value;
+	if (key.empty() || value.empty())
+		throw std::runtime_error("Cannot parse input line: " + line);
+	apply_key_value(config, current_test, key, value);
+}
+
+void parse_test_line(InputConfig& config, const std::string& line) {
+	std::istringstream in(line);
+	std::string marker;
+	std::string name;
+	in >> marker >> name;
+	if (name.empty())
+		throw std::runtime_error("Expected test name after '" + marker + "'.");
+
+	TestRunConfig& run = append_test(config, name);
+	std::string token;
+	while (in >> token) {
+		const size_t eq = token.find('=');
+		if (eq == std::string::npos)
+			throw std::runtime_error("Expected key=value in test line: " + token);
+		apply_override(run, token.substr(0, eq), token.substr(eq + 1));
+	}
+}
+
+InputConfig read_input(const std::filesystem::path& file_name) {
+	InputConfig config;
+	std::ifstream in(file_name);
+	if (!in.is_open()) {
+		std::cout << file_name.string()
+			<< " was not found; all tests will run with default parameters.\n";
+		add_test_name(config, "all");
+		return config;
+	}
+
+	TestRunConfig* current_test = nullptr;
+	std::string line;
+	size_t line_number = 0;
+	while (std::getline(in, line)) {
+		++line_number;
+		const size_t comment = line.find('#');
+		if (comment != std::string::npos)
+			line = line.substr(0, comment);
+
+		line = trim(line);
+		if (line.empty())
+			continue;
+
+		try {
+			if (line.front() == '[' && line.back() == ']') {
+				const std::string name = trim(line.substr(1, line.size() - 2));
+				current_test = &append_test(config, name);
+			}
+			else if (starts_with(normalize(line), "test ")) {
+				parse_test_line(config, line);
+				current_test = nullptr;
+			}
+			else {
+				parse_assignment_line(config, current_test, line);
+			}
+		}
+		catch (const std::exception& e) {
+			throw std::runtime_error(file_name.string() + ":" +
+				std::to_string(line_number) + ": " + e.what());
+		}
+	}
+
+	if (config.tests.empty())
+		add_test_name(config, "all");
+
+	return config;
+}
+
+MeshSize apply_overrides(const TestCase& test, const TestRunConfig& run) {
+	MeshSize mesh = test.mesh;
+	if (run.bottom_x) mesh.bottom_x = *run.bottom_x;
+	if (run.bottom_y) mesh.bottom_y = *run.bottom_y;
+	if (run.top_x) mesh.top_x = *run.top_x;
+	if (run.top_y) mesh.top_y = *run.top_y;
+	if (run.lambda_nodes) mesh.lambda_nodes = *run.lambda_nodes;
+	return mesh;
+}
+
+void validate_mesh(const MeshSize& mesh, ContactMethod contact_method) {
+	if (mesh.bottom_x < 2 || mesh.bottom_y < 2 || mesh.top_x < 2 || mesh.top_y < 2)
+		throw std::runtime_error("Every body mesh must have at least two nodes in each direction.");
+	if (contact_method != ContactMethod::SlaveNodes && mesh.lambda_nodes < 2)
+		throw std::runtime_error("lambda_nodes must be at least 2 for uniform contact methods.");
+}
+
+void validate_geometry(CoordinateSystem coordinate_system, const TestCase& test) {
+	if (coordinate_system == CoordinateSystem::Axisymmetric &&
+		(test.bottom_a.x <= 0.0 || test.top_a.x <= 0.0)) {
+		throw std::runtime_error("Axisymmetric tests require positive radial coordinates.");
+	}
+}
+
+std::string run_label(
+	CoordinateSystem coordinate_system,
+	ContactMethod contact_method,
+	ContactSlaveBody slave_body,
+	const MeshSize& mesh) {
+
+	std::ostringstream label;
+	label << "coord_" << to_string(coordinate_system)
+		<< "__contact_" << to_string(contact_method);
+	if (contact_method == ContactMethod::SlaveNodes)
+		label << "__passive_" << to_string(slave_body);
+	label << "__bottom_" << mesh.bottom_x << "x" << mesh.bottom_y
+		<< "__top_" << mesh.top_x << "x" << mesh.top_y;
+	if (contact_method != ContactMethod::SlaveNodes)
+		label << "__lambda_" << mesh.lambda_nodes;
+	return label.str();
+}
+
+double point_component(const Point& p, char component) {
+	return (component == 'x' || component == 'r') ? p.x : p.y;
+}
+
+double derivative_x(
+	const FEM& mesh,
+	const std::vector<Point>& field,
+	size_t row,
+	size_t col,
+	char component) {
+
+	const size_t mx = mesh.xsize();
+	auto value = [&](size_t c) {
+		return point_component(field[row * mx + c], component);
+	};
+
+	if (mx == 2) {
+		const double dx = mesh[row * mx + 1].x - mesh[row * mx].x;
+		return (value(1) - value(0)) / dx;
+	}
+
+	if (col == 0) {
+		const double dx = mesh[row * mx + 1].x - mesh[row * mx].x;
+		return (-3.0 * value(0) + 4.0 * value(1) - value(2)) / (2.0 * dx);
+	}
+
+	if (col + 1 == mx) {
+		const double dx = mesh[row * mx + col].x - mesh[row * mx + col - 1].x;
+		return (3.0 * value(col) - 4.0 * value(col - 1) + value(col - 2)) / (2.0 * dx);
+	}
+
+	const double dx = mesh[row * mx + col + 1].x - mesh[row * mx + col - 1].x;
+	return (value(col + 1) - value(col - 1)) / dx;
+}
+
+double derivative_y(
+	const FEM& mesh,
+	const std::vector<Point>& field,
+	size_t row,
+	size_t col,
+	char component) {
+
+	const size_t mx = mesh.xsize();
+	const size_t ny = mesh.ysize();
+	auto value = [&](size_t r) {
+		return point_component(field[r * mx + col], component);
+	};
+
+	if (ny == 2) {
+		const double dy = mesh[mx + col].y - mesh[col].y;
+		return (value(1) - value(0)) / dy;
+	}
+
+	if (row == 0) {
+		const double dy = mesh[mx + col].y - mesh[col].y;
+		return (-3.0 * value(0) + 4.0 * value(1) - value(2)) / (2.0 * dy);
+	}
+
+	if (row + 1 == ny) {
+		const double dy = mesh[row * mx + col].y - mesh[(row - 1) * mx + col].y;
+		return (3.0 * value(row) - 4.0 * value(row - 1) + value(row - 2)) / (2.0 * dy);
+	}
+
+	const double dy = mesh[(row + 1) * mx + col].y - mesh[(row - 1) * mx + col].y;
+	return (value(row + 1) - value(row - 1)) / dy;
+}
+
+double axisymmetric_hoop_strain(
+	const FEM& mesh,
+	const std::vector<Point>& field,
+	size_t row,
+	size_t col) {
+
+	const size_t node_id = row * mesh.xsize() + col;
+	const double r = mesh[node_id].x;
+	if (std::fabs(r) < kTraceEps)
+		return derivative_x(mesh, field, row, col, 'r');
+	return field[node_id].x / r;
+}
+
+std::vector<double> side_axis_coordinates(const FSEM& body, char side, bool fem_nodes) {
+	const auto nodes = fem_nodes ? body.get_side_fem_nodes(side) : body.get_side_nodes(side);
+	std::vector<double> x(nodes.size());
+	for (size_t i = 0; i < nodes.size(); ++i)
+		x[i] = fem_nodes ? body.fem[nodes[i]].x : body[nodes[i]].x;
+	return x;
+}
+
+size_t find_trace_segment(const std::vector<double>& x_nodes, double x) {
+	if (x_nodes.size() < 2)
+		throw std::runtime_error("At least two trace nodes are required.");
+
+	if (x <= x_nodes.front() + kTraceEps)
+		return 0;
+	if (x >= x_nodes.back() - kTraceEps)
+		return x_nodes.size() - 2;
+
+	for (size_t i = 0; i + 1 < x_nodes.size(); ++i)
+		if (x >= x_nodes[i] - kTraceEps && x <= x_nodes[i + 1] + kTraceEps)
+			return i;
+
+	throw std::runtime_error("Trace interpolation point is outside the contact interval.");
+}
+
+double interpolate_trace_value(
+	const std::vector<double>& x_nodes,
+	const std::vector<double>& values,
+	double x) {
+
+	const size_t segment = find_trace_segment(x_nodes, x);
+	const double x_left = x_nodes[segment];
+	const double x_right = x_nodes[segment + 1];
+	const double value_left = values[segment];
+	const double value_right = values[segment + 1];
+
+	if (std::fabs(x_left - x_right) < kTraceEps)
+		return value_left;
+
+	const double t = (x - x_left) / (x_right - x_left);
+	return (1.0 - t) * value_left + t * value_right;
+}
+
+std::vector<double> contact_output_grid(
+	const std::vector<double>& bottom_x,
+	const std::vector<double>& top_x) {
+
+	const double contact_left = std::max(bottom_x.front(), top_x.front());
+	const double contact_right = std::min(bottom_x.back(), top_x.back());
+
+	std::vector<double> grid;
+	grid.reserve(bottom_x.size() + top_x.size() + 2);
+	grid.push_back(contact_left);
+	grid.push_back(contact_right);
+
+	auto append = [&](const std::vector<double>& nodes) {
+		for (double x : nodes)
+			if (x >= contact_left - kTraceEps && x <= contact_right + kTraceEps)
+				grid.push_back(x);
+	};
+
+	append(bottom_x);
+	append(top_x);
+
+	std::sort(grid.begin(), grid.end());
+	grid.erase(std::unique(grid.begin(), grid.end(),
+		[](double lhs, double rhs) { return std::fabs(lhs - rhs) < kTraceEps; }),
+		grid.end());
+	return grid;
+}
+
+std::vector<double> recover_side_normal_stress(
+	CoordinateSystem coordinate_system,
+	const FSEM& body,
+	const std::vector<Point>& field,
+	char side,
+	double E,
+	double nu) {
+
+	const auto side_nodes = body.get_side_fem_nodes(side);
+	const double lambda = E * nu / ((1.0 + nu) * (1.0 - 2.0 * nu));
+	const double mu = E / (2.0 * (1.0 + nu));
+	const size_t mx = body.fem.xsize();
+
+	std::vector<double> sigma(side_nodes.size(), 0.0);
+	for (size_t i = 0; i < side_nodes.size(); ++i) {
+		const size_t node_id = side_nodes[i];
+		const size_t row = node_id / mx;
+		const size_t col = node_id % mx;
+
+		if (coordinate_system == CoordinateSystem::Axisymmetric) {
+			const double dur_dr = derivative_x(body.fem, field, row, col, 'r');
+			const double duz_dz = derivative_y(body.fem, field, row, col, 'z');
+			const double hoop_strain = axisymmetric_hoop_strain(body.fem, field, row, col);
+			sigma[i] = lambda * (dur_dr + hoop_strain) + (lambda + 2.0 * mu) * duz_dz;
+		}
+		else {
+			const double dux_dx = derivative_x(body.fem, field, row, col, 'x');
+			const double duy_dy = derivative_y(body.fem, field, row, col, 'y');
+			sigma[i] = lambda * dux_dx + (lambda + 2.0 * mu) * duy_dy;
+		}
+	}
+
+	return sigma;
+}
+
+std::vector<double> recover_side_normal_displacement(
+	const FSEM& body,
+	const std::vector<Point>& field,
+	char side) {
+
+	const auto side_nodes = body.get_side_fem_nodes(side);
+	std::vector<double> values(side_nodes.size(), 0.0);
+	for (size_t i = 0; i < side_nodes.size(); ++i)
+		values[i] = field[side_nodes[i]].y;
+	return values;
+}
+
+void save_displacement_component(
+	const std::filesystem::path& file_name,
+	const FEM& mesh,
+	const std::vector<Point>& field,
+	char component) {
+
+	std::ofstream out(file_name);
+	out << std::setprecision(16);
+	for (size_t i = 0; i < std::min(mesh.psize(), field.size()); ++i)
+		out << mesh[i].x << " " << mesh[i].y << " "
+			<< point_component(field[i], component) << "\n";
+}
+
+void save_displacement_vector(
+	const std::filesystem::path& file_name,
+	const FEM& mesh,
+	const std::vector<Point>& field) {
+
+	std::ofstream out(file_name);
+	out << std::setprecision(16);
+	for (size_t i = 0; i < std::min(mesh.psize(), field.size()); ++i)
+		out << mesh[i].x << " " << mesh[i].y << " " << field[i].x << " " << field[i].y << "\n";
+}
+
+void save_contact_normal_stress(
+	const std::filesystem::path& file_name,
+	CoordinateSystem coordinate_system,
+	const FSEM& bottom,
+	const std::vector<Point>& bottom_field,
+	const FSEM& top,
+	const std::vector<Point>& top_field,
+	double E,
+	double nu) {
+
+	const std::vector<double> bottom_x = side_axis_coordinates(bottom, 'N', true);
+	const std::vector<double> top_x = side_axis_coordinates(top, 'S', true);
+	const std::vector<double> grid = contact_output_grid(bottom_x, top_x);
+	const std::vector<double> bottom_sigma =
+		recover_side_normal_stress(coordinate_system, bottom, bottom_field, 'N', E, nu);
+	const std::vector<double> top_sigma =
+		recover_side_normal_stress(coordinate_system, top, top_field, 'S', E, nu);
+
+	std::ofstream out(file_name);
+	out << std::setprecision(16);
+	for (double x : grid)
+		out << x << " "
+			<< interpolate_trace_value(bottom_x, bottom_sigma, x) << " "
+			<< interpolate_trace_value(top_x, top_sigma, x) << "\n";
+}
+
+void save_contact_normal_displacement(
+	const std::filesystem::path& file_name,
+	const FSEM& bottom,
+	const std::vector<Point>& bottom_field,
+	const FSEM& top,
+	const std::vector<Point>& top_field) {
+
+	const std::vector<double> bottom_x = side_axis_coordinates(bottom, 'N', true);
+	const std::vector<double> top_x = side_axis_coordinates(top, 'S', true);
+	const std::vector<double> grid = contact_output_grid(bottom_x, top_x);
+	const std::vector<double> bottom_u = recover_side_normal_displacement(bottom, bottom_field, 'N');
+	const std::vector<double> top_u = recover_side_normal_displacement(top, top_field, 'S');
+
+	std::ofstream out(file_name);
+	out << std::setprecision(16);
+	for (double x : grid) {
+		const double ub = interpolate_trace_value(bottom_x, bottom_u, x);
+		const double ut = interpolate_trace_value(top_x, top_u, x);
+		out << x << " " << ub << " " << ut << " " << (ub - ut) << "\n";
+	}
+}
+
+std::vector<double> lambda_nodes_for_output(
+	const FSEM& bottom,
+	const FSEM& top,
+	const ContactOptions& options) {
+
+	const std::vector<double> bottom_x = side_axis_coordinates(bottom, 'N', false);
+	const std::vector<double> top_x = side_axis_coordinates(top, 'S', false);
+	const double contact_left = std::max(bottom_x.front(), top_x.front());
+	const double contact_right = std::min(bottom_x.back(), top_x.back());
+
+	ContactDiscretization discretization;
+	if (options.method == ContactMethod::SlaveNodes) {
+		discretization = contact_slave_nodes::build(
+			bottom_x, top_x, contact_left, contact_right, options.slave_body);
+	}
+	else if (options.method == ContactMethod::UniformLambdaPartition) {
+		discretization = contact_uniform_lambda_partition::build(
+			bottom_x, top_x, contact_left, contact_right, options.lambda_node_count);
+	}
+	else {
+		discretization = contact_uniform_union_partition::build(
+			bottom_x, top_x, contact_left, contact_right, options.lambda_node_count);
+	}
+	return discretization.lambda_nodes;
+}
+
+void save_lagrange_multipliers(
+	const std::filesystem::path& file_name,
+	const std::vector<double>& solution,
+	size_t start,
+	const std::vector<double>& lambda_nodes) {
+
+	std::ofstream out(file_name);
+	out << std::setprecision(16);
+	const size_t count = std::min(lambda_nodes.size(), solution.size() - start);
+	for (size_t i = 0; i < count; ++i)
+		out << lambda_nodes[i] << " " << solution[start + i] << "\n";
+}
+
+double relative_error(
+	const FEM& mesh,
+	const std::vector<Point>& field,
+	const vec_function& exact) {
+
+	double numerator = 0.0;
+	double denominator = 0.0;
+	for (size_t i = 0; i < std::min(mesh.psize(), field.size()); ++i) {
+		const Point u = exact(mesh[i]);
+		const Point diff = field[i] - u;
+		numerator += diff.x * diff.x + diff.y * diff.y;
+		denominator += u.x * u.x + u.y * u.y;
+	}
+
+	return denominator > 1e-30 ? std::sqrt(numerator / denominator) : std::sqrt(numerator);
+}
+
+void save_parameters(
+	const std::filesystem::path& file_name,
+	const std::string& test_name,
+	const MeshSize& mesh,
+	const ContactOptions& options,
+	double E,
+	double nu) {
+
+	std::ofstream out(file_name);
+	out << "test=" << test_name << "\n";
+	out << "coordinate=" << to_string(options.coordinate_system) << "\n";
+	out << "contact=" << to_string(options.method) << "\n";
+	out << "passive_body=" << to_string(options.slave_body) << "\n";
+	out << "bottom_mesh=" << mesh.bottom_x << "x" << mesh.bottom_y << "\n";
+	out << "top_mesh=" << mesh.top_x << "x" << mesh.top_y << "\n";
+	out << "lambda_nodes=" << mesh.lambda_nodes << "\n";
+	out << "E=" << std::setprecision(16) << E << "\n";
+	out << "nu=" << nu << "\n";
+}
+
+void apply_dirichlet(FSEM& body, const std::array<char, 3>& sides, const vec_function& exact) {
+	for (char side : sides)
+		body.set_bc1(side, exact);
+}
+
+void run_test(const TestRunConfig& run, const InputConfig& config) {
+	const auto it = TESTS.find(run.name);
+	if (it == TESTS.end()) {
+		std::cerr << "Unknown test '" << run.name << "'. It is skipped.\n";
+		return;
+	}
+
+	const TestCase& test = it->second;
+	const MeshSize mesh = apply_overrides(test, run);
+	validate_mesh(mesh, config.contact_method);
+	validate_geometry(config.coordinate_system, test);
+
+	const vec_function exact = test.exact_solution(test.nu, test.E);
+	FSEM bottom(test.E, test.nu, test.bottom_a, test.bottom_b,
+		mesh.bottom_x, mesh.bottom_y, 1, 1, config.coordinate_system);
+	FSEM top(test.E, test.nu, test.top_a, test.top_b,
+		mesh.top_x, mesh.top_y, 1, 1, config.coordinate_system);
+
+	bottom.construct_basis();
+	top.construct_basis();
+	apply_dirichlet(bottom, test.bottom_dirichlet_sides, exact);
+	apply_dirichlet(top, test.top_dirichlet_sides, exact);
+
+	ContactOptions contact_options;
+	contact_options.coordinate_system = config.coordinate_system;
+	contact_options.method = config.contact_method;
+	contact_options.slave_body = config.slave_body;
+	contact_options.lambda_node_count = mesh.lambda_nodes;
+
+	const std::vector<double> solution = solve_mortar_contact(
+		bottom, top, bottom.get_f(), top.get_f(), contact_options);
+
+	const size_t n_bottom = bottom.get_K().size();
+	const size_t n_top = top.get_K().size();
+	const std::vector<Point> bottom_field = bottom.find_answer(solution);
+	const std::vector<Point> top_field = top.find_answer(solution, n_bottom);
+
+	const std::filesystem::path output_dir =
+		std::filesystem::path("res") / run.name /
+		run_label(config.coordinate_system, config.contact_method, config.slave_body, mesh);
+	std::filesystem::create_directories(output_dir);
+
+	const char first_component = config.coordinate_system == CoordinateSystem::Axisymmetric ? 'r' : 'x';
+	const char second_component = config.coordinate_system == CoordinateSystem::Axisymmetric ? 'z' : 'y';
+
+	save_displacement_vector(output_dir / "bottom_displacement.txt", bottom.fem, bottom_field);
+	save_displacement_vector(output_dir / "top_displacement.txt", top.fem, top_field);
+	save_displacement_component(output_dir / ("bottom_displacement_" + std::string(1, first_component) + ".txt"),
+		bottom.fem, bottom_field, first_component);
+	save_displacement_component(output_dir / ("bottom_displacement_" + std::string(1, second_component) + ".txt"),
+		bottom.fem, bottom_field, second_component);
+	save_displacement_component(output_dir / ("top_displacement_" + std::string(1, first_component) + ".txt"),
+		top.fem, top_field, first_component);
+	save_displacement_component(output_dir / ("top_displacement_" + std::string(1, second_component) + ".txt"),
+		top.fem, top_field, second_component);
+
+	save_contact_normal_stress(output_dir / "contact_normal_stress.txt",
+		config.coordinate_system, bottom, bottom_field, top, top_field, test.E, test.nu);
+	save_contact_normal_displacement(output_dir / "contact_normal_displacement.txt",
+		bottom, bottom_field, top, top_field);
+	save_lagrange_multipliers(output_dir / "lagrange_multipliers.txt",
+		solution, n_bottom + n_top, lambda_nodes_for_output(bottom, top, contact_options));
+
+	const double bottom_error = relative_error(bottom.fem, bottom_field, exact);
+	const double top_error = relative_error(top.fem, top_field, exact);
+	{
+		std::ofstream out(output_dir / "error.txt");
+		out << std::setprecision(16)
+			<< "bottom_relative_error=" << bottom_error << "\n"
+			<< "top_relative_error=" << top_error << "\n";
+	}
+	save_parameters(output_dir / "parameters.txt", run.name, mesh, contact_options, test.E, test.nu);
+
+	std::cout << "Saved " << run.name << " -> " << output_dir.string()
+		<< " (lambda=" << (solution.size() - n_bottom - n_top)
+		<< ", error bottom=" << bottom_error
+		<< ", top=" << top_error << ")\n";
+}
+
+} // namespace
+
+int main(int argc, char* argv[]) {
+	try {
+		if (argc > 2) {
+			std::cerr << "Usage: mkse-elasticity [input-file]\n";
+			return 2;
+		}
+
+		const std::filesystem::path input_file = argc == 2 ? argv[1] : "input.txt";
+		const InputConfig config = read_input(input_file);
+		std::cout << "Coordinate system: " << to_string(config.coordinate_system) << "\n";
+		std::cout << "Contact method: " << to_string(config.contact_method) << "\n";
+		if (config.contact_method == ContactMethod::SlaveNodes)
+			std::cout << "Passive body: " << to_string(config.slave_body) << "\n";
+
+		for (const TestRunConfig& run : config.tests) {
+			try {
+				run_test(run, config);
+			}
+			catch (const std::exception& e) {
+				std::cerr << "Test '" << run.name << "' failed: " << e.what() << "\n";
+			}
+		}
+	}
+	catch (const std::exception& e) {
+		std::cerr << "Fatal error: " << e.what() << "\n";
+		return 1;
+	}
+
+	return 0;
+}