diff --git a/.gitattributes b/.gitattributes new file mode 100644 index 0000000..552b608 --- /dev/null +++ b/.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 diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..8bda51e --- /dev/null +++ b/.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 diff --git a/CMakeLists.txt b/CMakeLists.txt new file mode 100644 index 0000000..790a6b1 --- /dev/null +++ b/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() diff --git a/include/ContactSlaveNodes.h b/include/ContactSlaveNodes.h new file mode 100644 index 0000000..57559e2 --- /dev/null +++ b/include/ContactSlaveNodes.h @@ -0,0 +1,91 @@ +#pragma once + +// Contact discretization based on the passive body's trace nodes. + +#include +#include +#include + +#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 collect_overlap_nodes( + const std::vector& bottom_x, + const std::vector& top_x, + double contact_left, + double contact_right) { + + std::vector 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& 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 collect_active_lambda_nodes( + const std::vector& passive_x, + double contact_left, + double contact_right) { + + std::vector 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& bottom_x, + const std::vector& top_x, + double contact_left, + double contact_right, + ContactSlaveBody passive_body) { + + const std::vector& 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 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 diff --git a/include/ContactUniformLambdaPartition.h b/include/ContactUniformLambdaPartition.h new file mode 100644 index 0000000..3247c01 --- /dev/null +++ b/include/ContactUniformLambdaPartition.h @@ -0,0 +1,81 @@ +#pragma once + +// Contact discretization on an independent uniform multiplier grid. + +#include +#include + +#include "Elasticity.h" + +namespace contact_uniform_lambda_partition { + +constexpr double kContactEps = 1e-12; + +inline size_t count_contact_nodes( + const std::vector& 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& bottom_x, + const std::vector& top_x, + double contact_left, + double contact_right) { + + return std::max( + 2, + std::max( + count_contact_nodes(bottom_x, contact_left, contact_right), + count_contact_nodes(top_x, contact_left, contact_right))); +} + +inline std::vector build_uniform_lambda_nodes( + double contact_left, + double contact_right, + size_t lambda_node_count) { + + std::vector 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& bottom_x, + const std::vector& 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(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 diff --git a/include/ContactUniformUnionPartition.h b/include/ContactUniformUnionPartition.h new file mode 100644 index 0000000..03ebf78 --- /dev/null +++ b/include/ContactUniformUnionPartition.h @@ -0,0 +1,73 @@ +#pragma once + +// Contact discretization on a uniform refinement of the combined traces. + +#include +#include + +#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 build_mortar_elements( + const std::vector& bottom_x, + const std::vector& top_x, + const std::vector& lambda_nodes, + double contact_left, + double contact_right) { + + std::vector 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& 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 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& bottom_x, + const std::vector& 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 diff --git a/include/Elasticity.h b/include/Elasticity.h new file mode 100644 index 0000000..441807d --- /dev/null +++ b/include/Elasticity.h @@ -0,0 +1,265 @@ +#pragma once + +// Core data structures and solvers for the finite superelement model. + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +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 operator-(const std::vector& a, + const std::vector& b); + +using function = std::function; +using vec_function = std::function; +using lambda_func = std::function; + +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 lambda_nodes; + std::vector mortar_elements; +}; + +class Matrix { +private: + std::vector> matrix; + +public: + Matrix(size_t n, double a = 0); + Matrix(size_t n, size_t m, double a = 0); + Matrix(std::vector> m) : matrix(m) {}; + Matrix(Matrix* M) { matrix = M->matrix; }; + + std::vector& operator[](size_t i) { return matrix[i]; }; + const std::vector& 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 dot(const std::vector& v) const; + void print() const; + + static Matrix eye(size_t n, double a = 1); +}; + +std::vector solveGaussFullPivot( + const Matrix& A, + const std::vector& b, + double eps = 1e-12 +); + +class FEM { + size_t mx, ny; + Point left_down; + Point right_up; + std::vector points; + std::vector triangles; + double triangle_area; + std::vector f; + std::vector u; + CoordinateSystem coordinate_system; + + std::vector 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 solve(); + + void clear_AFu(); + + void bc2_side(lambda_func j, size_t start, size_t finish, double len, + int side, const std::vector& g, + std::vector& p_vec); + + void calculate_bc2(const std::vector& pos, + const std::vector& g, std::vector& p_vec); + + std::pair> get_AF(); + void set_AF(const Matrix& A_new, const std::vector& F_new); +}; + +Matrix operator*(double a, Matrix m); + +std::tuple LU_decomposition(const Matrix& m); + +std::vector solveLU(const Matrix& L, const Matrix& U, + const std::vector& 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 nodes; + size_t coef_x; + size_t coef_y; + + std::vector> basis; + std::vector basis_coefficients; + + Matrix K; + std::vector 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& get_f() const { return f; } + const std::vector>& get_basis() const { return basis; } + + void construct_f_bc2(const std::vector& pos, + const std::vector& g); + + std::vector find_answer(); + std::vector find_answer(const std::vector& 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& coefs_Dirichle, int& dir_id, const int i); + + void set_bc2(const std::vector& pos, + const std::vector& g); + + std::vector get_side_nodes(char side) const; + std::vector 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> get_known_dofs() const; +}; + +double mortar_shape_func(size_t i, const std::vector& s, double cur); + +std::vector solve_mortar_contact( + FSEM& bottom_body, + FSEM& top_body, + const std::vector& rhs_bottom, + const std::vector& rhs_top, + const ContactOptions& options); + +std::vector solve_mortar_contact( + FSEM& bottom_body, + FSEM& top_body, + const std::vector& rhs_bottom, + const std::vector& rhs_top, + size_t lambda_node_count = 0); + +std::vector solveWithLU(const Matrix& A, + const std::vector& b, + double eps = 1e-15); diff --git a/include/FEMAxisymmetric.h b/include/FEMAxisymmetric.h new file mode 100644 index 0000000..e1e8e0f --- /dev/null +++ b/include/FEMAxisymmetric.h @@ -0,0 +1,115 @@ +#pragma once + +// Axisymmetric finite-element integration kernels. + +#include +#include + +#include "Elasticity.h" + +namespace fem_axisymmetric { + +constexpr double kTwoPi = 6.28318530717958647692; + +struct TriangleQuadraturePoint { + double weight; + std::array phi; +}; + +inline const std::array 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& F, + const std::vector& points, + const std::vector& 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 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 diff --git a/include/FEMCartesian.h b/include/FEMCartesian.h new file mode 100644 index 0000000..e185d03 --- /dev/null +++ b/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& F, + const std::vector& points, + const std::vector& 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 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 diff --git a/include/TestCases.h b/include/TestCases.h new file mode 100644 index 0000000..a1a2205 --- /dev/null +++ b/include/TestCases.h @@ -0,0 +1,114 @@ +#pragma once + +// Manufactured solutions used to validate the numerical model. + +#include +#include +#include +#include +#include + +#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 exact_solution; + double E; + double nu; + Point bottom_a; + Point bottom_b; + Point top_a; + Point top_b; + MeshSize mesh; + std::array bottom_dirichlet_sides; + std::array top_dirichlet_sides; +}; + +const std::map 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{ 'W', 'E', 'S' }, + std::array{ '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{ 'W', 'E', 'S' }, + std::array{ '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{ 'W', 'E', 'S' }, + std::array{ '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{ 'W', 'E', 'S' }, + std::array{ 'W', 'N', 'E' } + } + } +}; diff --git a/input.txt b/input.txt new file mode 100644 index 0000000..4f7d37e --- /dev/null +++ b/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 \ No newline at end of file diff --git a/mkse-elasticity.sln b/mkse-elasticity.sln new file mode 100644 index 0000000..5278706 --- /dev/null +++ b/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 diff --git a/mkse-elasticity.vcxproj b/mkse-elasticity.vcxproj new file mode 100644 index 0000000..09e6c29 --- /dev/null +++ b/mkse-elasticity.vcxproj @@ -0,0 +1,165 @@ + + + + + Debug + Win32 + + + Release + Win32 + + + Debug + x64 + + + Release + x64 + + + + 16.0 + Win32Proj + {9e82b28a-0575-4147-bc8d-b9e212c23900} + mkse_elasticity + 10.0 + + + + Application + true + v143 + Unicode + + + Application + false + v143 + true + Unicode + + + Application + true + v143 + Unicode + + + Application + false + v143 + true + Unicode + + + + + + + + + + + + + + + + + + + + + $(ProjectDir)build\msbuild\bin\$(Platform)\$(Configuration)\ + $(ProjectDir)build\msbuild\obj\$(Platform)\$(Configuration)\ + mkse-elasticity + + + + $(ProjectDir)include;%(AdditionalIncludeDirectories) + + + + + Level4 + true + WIN32;_DEBUG;_CONSOLE;%(PreprocessorDefinitions) + true + stdcpp20 + + + Console + true + + + + + Level4 + true + true + true + WIN32;NDEBUG;_CONSOLE;%(PreprocessorDefinitions) + true + stdcpp20 + + + Console + true + true + true + + + + + Level4 + true + _DEBUG;_CONSOLE;%(PreprocessorDefinitions) + true + stdcpp20 + + + Console + true + + + + + Level4 + true + true + true + NDEBUG;_CONSOLE;%(PreprocessorDefinitions) + true + stdcpp20 + + + Console + true + true + true + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/mkse-elasticity.vcxproj.filters b/mkse-elasticity.vcxproj.filters new file mode 100644 index 0000000..1d9494a --- /dev/null +++ b/mkse-elasticity.vcxproj.filters @@ -0,0 +1,54 @@ + + + + + {4FC737F1-C7A5-4376-A066-2A32D752A2FF} + cpp;c;cc;cxx;c++;cppm;ixx;def;odl;idl;hpj;bat;asm;asmx + + + {93995380-89BD-4b04-88EB-625FBE52EBFB} + h;hh;hpp;hxx;h++;hm;inl;inc;ipp;xsd + + + {67DA6AB6-F800-4c08-8B7A-83BB121AAD01} + rc;ico;cur;bmp;dlg;rc2;rct;bin;rgs;gif;jpg;jpeg;jpe;resx;tiff;tif;png;wav;mfcribbon-ms + + + + + Source Files + + + Source Files + + + Source Files + + + Source Files + + + + + Header Files + + + Header Files + + + Header Files + + + Header Files + + + Header Files + + + Header Files + + + Header Files + + + diff --git a/src/FEM.cpp b/src/FEM.cpp new file mode 100644 index 0000000..e3447a2 --- /dev/null +++ b/src/FEM.cpp @@ -0,0 +1,274 @@ +#include +#include +#include +#include +#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 phi; + }; + + const std::array 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(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 FEM::solve() { + auto [L, U] = LU_decomposition(A); + std::vector dofs = solveLU(L, U, F); + + std::vector 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> FEM::get_AF() { + return { A, F }; +} + +void FEM::set_AF(const Matrix& A_new, const std::vector& 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& g, + std::vector& 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& pos, + const std::vector& g, std::vector& 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); + +} diff --git a/src/FSEM.cpp b/src/FSEM.cpp new file mode 100644 index 0000000..17619a4 --- /dev/null +++ b/src/FSEM.cpp @@ -0,0 +1,1097 @@ +#include "Elasticity.h" +#include "ContactSlaveNodes.h" +#include "ContactUniformLambdaPartition.h" +#include "ContactUniformUnionPartition.h" +#include +#include + +Point zero(const Point& p) { + (void)p; + return { 0, 0 }; +} + +std::string to_string(CoordinateSystem coordinate_system) { + return coordinate_system == CoordinateSystem::Cartesian + ? "cartesian" + : "axisymmetric"; +} + +std::string to_string(ContactMethod contact_method) { + if (contact_method == ContactMethod::SlaveNodes) + return "slave-nodes"; + if (contact_method == ContactMethod::UniformLambdaPartition) + return "uniform-lambda"; + return "uniform-union"; +} + +std::string to_string(ContactSlaveBody slave_body) { + return slave_body == ContactSlaveBody::Bottom ? "bottom" : "top"; +} + + +//FSEM::FSEM(double E, double nu, const Point& a, const Point& b, +// size_t n_x, size_t n_y, int coef_val_x, int coef_val_y): +// E(E), nu(nu), a(a), b(b), n_side_x(n_x - 1), +// n_side_y(n_y - 1), coef_x(coef_val_x), coef_y(coef_val_y), +// fem(a, b, (n_y - 1)* coef_y + 1, (n_x - 1)* coef_x + 1) { +FSEM::FSEM(double E, double nu, const Point& a, const Point& b, + size_t n_x, size_t n_y, int coef_val_x, int coef_val_y, + CoordinateSystem coordinate_system) : + a(a), b(b), E(E), nu(nu), n_side_x(n_x - 1), + n_side_y(n_y - 1), coef_x(static_cast(coef_val_x)), + coef_y(static_cast(coef_val_y)), + K(2 * (coef_val_x * (n_x - 1) + 1) * (coef_val_y * (n_y - 1) + 1)), + f(2 * (coef_val_x * (n_x - 1) + 1) * (coef_val_y * (n_y - 1) + 1)), + coordinate_system(coordinate_system), + fem(a, b, (n_y - 1) * coef_val_y + 1, (n_x - 1) * coef_val_x + 1, + coordinate_system) + { + + double h_x = (b.x - a.x) / n_side_x, h_y = (b.y - a.y) / n_side_y; // шаг + + const int n_nodes = static_cast(2 * (n_side_x + n_side_y)); + nodes.resize(n_nodes); + + basis = std::vector>(2 * n_nodes, + std::vector((coef_x * n_side_x + 1) * (coef_y * n_side_y + 1))); + + basis_coefficients = std::vector(2 * n_nodes, { NAN, NAN }); + + for (int i = 0; i < n_side_y; ++i) { + nodes[i] = { a.x, a.y + i * h_y }; + nodes[n_side_x + n_side_y + i] = { b.x, b.y - i * h_y }; + } + for (int i = 0; i < n_side_x; ++i) { + nodes[n_side_y + i] = { a.x + i * h_x, b.y }; + nodes[n_side_x + 2 * n_side_y + i] = { b.x - i * h_x, a.y }; + } + +} + +Point& FSEM::get_node(size_t n) { + return nodes[n]; +} + +const Point& FSEM::get_node(size_t n) const { + return nodes[n]; +} + +void FSEM::print_nodes() const { + for (size_t i = 0; i < nodes.size(); ++i) { + std::cout << (*this)[i] << "\t"; + } + std::cout << '\n'; +} + +void FSEM::construct_basis() { + // количество узлов сетки для МКЭ + /*int n_x = n_side_x * coef_x + 1, n_y = n_side_y * coef_y + 1; + double dx = (b.x - a.x) / (n_x - 1), dy = (b.y - a.y) / (n_y - 1);*/ + double dx = (b.x - a.x) / n_side_x, dy = (b.y - a.y) / n_side_y; // шаг + + //******************ЛЕВЫЙ НИЖНИЙ УГОЛ****************** + + // находим базисную функцию для, равную 1 в a по х компоненте + fem.set_boundaries('W', + [&](const Point& p) + { + double neighbor_y = a.y + dy; + if (p.y < neighbor_y) + return Point{ (neighbor_y - p.y) / dy, 0 }; + + return Point{ 0, 0 }; + }); + fem.set_boundaries('N', zero); + fem.set_boundaries('E', zero); + fem.set_boundaries('S', + [&](const Point& p) + { + double neighbor_x = a.x + dx; + if (p.x < neighbor_x) + return Point{ (neighbor_x - p.x) / dx, 0 }; + + return Point{ 0, 0 }; + }); + + fem.construct_AF(E, nu, zero); + auto [A, F] = fem.get_AF(); + fem.apply_boundaries(); + + basis[0] = fem.solve(); + + // находим базисную функцию для, равную 1 в a по y компоненте + fem.set_boundaries('W', + [&](const Point& p) + { + double neighbor_y = a.y + dy; + if (p.y < neighbor_y) + return Point{ 0, (neighbor_y - p.y) / dy}; + + return Point{ 0, 0 }; + }); + fem.set_boundaries('N', zero); + fem.set_boundaries('E', zero); + fem.set_boundaries('S', + [&](const Point& p) + { + double neighbor_x = a.x + dx; + if (p.x < neighbor_x) + return Point{ 0, (neighbor_x - p.x) / dx }; + + return Point{ 0, 0 }; + }); + fem.set_AF(A, F); + fem.apply_boundaries(); + + basis[1] = fem.solve(); + + //******************ЛЕВЫЙ ВЕРХНИЙ УГОЛ****************** + + // находим базисную функцию для, равную 1 в {a.x, b.y} по х компоненте + fem.set_boundaries('W', + [&](const Point& p) + { + double neighbor_y = b.y - dy; + if (p.y > neighbor_y) + return Point{ (p.y - neighbor_y) / dy, 0 }; + + return Point{ 0, 0 }; + }); + fem.set_boundaries('N', + [&](const Point& p) + { + double neighbor_x = a.x + dx; + if (p.x < neighbor_x) + return Point{ (neighbor_x - p.x) / dx, 0 }; + + return Point{ 0, 0 }; + }); + fem.set_boundaries('E', zero); + fem.set_boundaries('S', zero); + + fem.set_AF(A, F); + fem.apply_boundaries(); + basis[2 * n_side_y] = fem.solve(); + + // находим базисную функцию для, равную 1 в {a.x, b.y} по y компоненте + fem.set_boundaries('W', + [&](const Point& p) + { + double neighbor_y = b.y - dy; + if (p.y > neighbor_y) + return Point{ 0, (p.y - neighbor_y) / dy}; + + return Point{ 0, 0 }; + }); + fem.set_boundaries('N', + [&](const Point& p) + { + double neighbor_x = a.x + dx; + if (p.x < neighbor_x) + return Point{ 0, (neighbor_x - p.x) / dx}; + + return Point{ 0, 0 }; + }); + fem.set_boundaries('E', zero); + fem.set_boundaries('S', zero); + + fem.set_AF(A, F); + fem.apply_boundaries(); + basis[2 * n_side_y + 1] = fem.solve(); + + //******************ПРАВЫЙ ВЕРХНИЙ УГОЛ****************** + + // находим базисную функцию для, равную 1 в b по х компоненте + fem.set_boundaries('W', zero ); + fem.set_boundaries('N', + [&](const Point& p) + { + double neighbor_x = b.x - dx; + if (p.x > neighbor_x) + return Point{ (p.x - neighbor_x) / dx, 0 }; + + return Point{ 0, 0 }; + }); + fem.set_boundaries('E', + [&](const Point& p) + { + double neighbor_y = b.y - dy; + if (p.y > neighbor_y) + return Point{ (p.y - neighbor_y) / dy, 0 }; + + return Point{ 0, 0 }; + }); + fem.set_boundaries('S', zero); + + fem.set_AF(A, F); + fem.apply_boundaries(); + basis[2 * (n_side_x + n_side_y)] = fem.solve(); + + // находим базисную функцию для, равную 1 в b по y компоненте + fem.set_boundaries('W', zero); + fem.set_boundaries('N', + [&](const Point& p) + { + double neighbor_x = b.x - dx; + if (p.x > neighbor_x) + return Point{ 0, (p.x - neighbor_x) / dx}; + + return Point{ 0, 0 }; + }); + fem.set_boundaries('E', + [&](const Point& p) + { + double neighbor_y = b.y - dy; + if (p.y > neighbor_y) + return Point{ 0, (p.y - neighbor_y) / dy}; + + return Point{ 0, 0 }; + }); + fem.set_boundaries('S', zero); + + fem.set_AF(A, F); + fem.apply_boundaries(); + basis[2 * (n_side_x + n_side_y) + 1] = fem.solve(); + + //******************ПРАВЫЙ НИЖНИЙ УГОЛ****************** + + // находим базисную функцию для, равную 1 в {b.x, a.y} по х компоненте + fem.set_boundaries('W', zero); + fem.set_boundaries('N', zero); + fem.set_boundaries('E', + [&](const Point& p) + { + double neighbor_y = a.y + dy; + if (p.y < neighbor_y) + return Point{ (neighbor_y - p.y) / dy, 0 }; + + return Point{ 0, 0 }; + }); + fem.set_boundaries('S', + [&](const Point& p) + { + double neighbor_x = b.x - dx; + if (p.x > neighbor_x) + return Point{ (p.x - neighbor_x) / dx, 0 }; + + return Point{ 0, 0 }; + }); + + fem.set_AF(A, F); + fem.apply_boundaries(); + basis[2 * (n_side_x + 2 * n_side_y)] = fem.solve(); + + // находим базисную функцию для, равную 1 в {b.x, a.y} по y компоненте + fem.set_boundaries('W', zero); + fem.set_boundaries('N', zero); + fem.set_boundaries('E', + [&](const Point& p) + { + double neighbor_y = a.y + dy; + if (p.y < neighbor_y) + return Point{ 0, (neighbor_y - p.y) / dy }; + + return Point{ 0, 0 }; + }); + fem.set_boundaries('S', + [&](const Point& p) + { + double neighbor_x = b.x - dx; + if (p.x > neighbor_x) + return Point{ 0, (p.x - neighbor_x) / dx}; + + return Point{ 0, 0 }; + }); + + fem.set_AF(A, F); + fem.apply_boundaries(); + basis[2 * (n_side_x + 2 * n_side_y) + 1] = fem.solve(); + + //******************ВЕРТИКАЛЬНАЯ СТОРОНА****************** + + for (int i = 1; i < n_side_y; ++i) { + + //--------------НА ЛЕВОЙ-------------- + + // находим базисную функцию для, равную 1 в i по х компоненте + fem.set_boundaries('W', + [&](const Point& p) + { + double current_y = a.y + i * dy; + double neighbor_up_y = current_y + dy, + neighbor_down_y = current_y - dy; + + if (neighbor_down_y < p.y && p.y <= current_y ) + return Point{ (p.y - neighbor_down_y) / dy, 0 }; + + if ( current_y < p.y && p.y < neighbor_up_y ) + return Point{ (neighbor_up_y - p.y) / dy, 0 }; + + return Point{ 0, 0 }; + }); + fem.set_boundaries('N', zero); + fem.set_boundaries('E', zero); + fem.set_boundaries('S', zero); + + fem.set_AF(A, F); + fem.apply_boundaries(); + basis[2 * i] = fem.solve(); + + // находим базисную функцию для, равную 1 в i по y компоненте + fem.set_boundaries('W', + [&](const Point& p) + { + double current_y = a.y + i * dy; + double neighbor_up_y = current_y + dy, + neighbor_down_y = current_y - dy; + + if (neighbor_down_y < p.y && p.y <= current_y) + return Point{ 0, (p.y - neighbor_down_y) / dy}; + + if (current_y < p.y && p.y < neighbor_up_y) + return Point{ 0, (neighbor_up_y - p.y) / dy}; + + return Point{ 0, 0 }; + }); + fem.set_boundaries('N', zero); + fem.set_boundaries('E', zero); + fem.set_boundaries('S', zero); + + fem.set_AF(A, F); + fem.apply_boundaries(); + basis[2 * i + 1] = fem.solve(); + + //--------------НА ПРАВОЙ-------------- + + // находим базисную функцию для, равную 1 в (n_side_x + n_side_y + i) по х компоненте + fem.set_boundaries('W', zero); + fem.set_boundaries('N', zero); + fem.set_boundaries('E', + [&](const Point& p) + { + double current_y = b.y - i * dy; + double neighbor_up_y = current_y + dy, + neighbor_down_y = current_y - dy; + + if (neighbor_down_y < p.y && p.y <= current_y) + return Point{ (p.y - neighbor_down_y) / dy, 0 }; + + if (current_y < p.y && p.y < neighbor_up_y) + return Point{ (neighbor_up_y - p.y) / dy, 0 }; + + return Point{ 0, 0 }; + }); + fem.set_boundaries('S', zero); + + fem.set_AF(A, F); + fem.apply_boundaries(); + basis[2 * (n_side_x + n_side_y + i)] = fem.solve(); + + // находим базисную функцию для, равную 1 в i по y компоненте + fem.set_boundaries('W', zero); + fem.set_boundaries('N', zero); + fem.set_boundaries('E', + [&](const Point& p) + { + double current_y = b.y - i * dy; + double neighbor_up_y = current_y + dy, + neighbor_down_y = current_y - dy; + + if (neighbor_down_y < p.y && p.y <= current_y) + return Point{ 0, (p.y - neighbor_down_y) / dy}; + + if (current_y < p.y && p.y < neighbor_up_y) + return Point{ 0, (neighbor_up_y - p.y) / dy }; + + return Point{ 0, 0 }; + }); + fem.set_boundaries('S', zero); + + fem.set_AF(A, F); + fem.apply_boundaries(); + basis[2 * (n_side_x + n_side_y + i) + 1] = fem.solve(); + } + + //******************ГОРИЗОНТАЛЬНАЯ СТОРОНА****************** + + for (int i = 1; i < n_side_x; ++i) { + + //--------------НА ВЕРХНЕЙ-------------- + + // находим базисную функцию для, равную 1 в n_side_x + i по х компоненте + fem.set_boundaries('W', zero); + fem.set_boundaries('N', + [&](const Point& p) + { + double current_x = a.x + i * dx; + double neighbor_right_x = current_x + dx, + neighbor_left_x = current_x - dx; + + if (neighbor_left_x < p.x && p.x <= current_x) + return Point{ (p.x - neighbor_left_x) / dx, 0 }; + + if (current_x < p.x && p.x < neighbor_right_x) + return Point{ (neighbor_right_x - p.x) / dx, 0 }; + + return Point{ 0, 0 }; + }); + fem.set_boundaries('E', zero); + fem.set_boundaries('S', zero); + + fem.set_AF(A, F); + fem.apply_boundaries(); + basis[2 * (n_side_y + i)] = fem.solve(); + + // находим базисную функцию для, равную 1 в n_side_x + i по y компоненте + fem.set_boundaries('W', zero); + fem.set_boundaries('N', + [&](const Point& p) + { + double current_x = a.x + i * dx; + double neighbor_right_x = current_x + dx, + neighbor_left_x = current_x - dx; + + if (neighbor_left_x < p.x && p.x <= current_x) + return Point{ 0, (p.x - neighbor_left_x) / dx}; + + if (current_x < p.x && p.x < neighbor_right_x) + return Point{ 0, (neighbor_right_x - p.x) / dx }; + + return Point{ 0, 0 }; + }); + fem.set_boundaries('E', zero); + fem.set_boundaries('S', zero); + + fem.set_AF(A, F); + fem.apply_boundaries(); + basis[2 * (n_side_y + i) + 1] = fem.solve(); + + //--------------НА НИЖНЕЙ-------------- + + // находим базисную функцию для, равную 1 в (2 * n_side_y + n_side_x + i) по х компоненте + fem.set_boundaries('W', zero); + fem.set_boundaries('N', zero); + fem.set_boundaries('E', zero); + fem.set_boundaries('S', + [&](const Point& p) + { + double current_x = b.x - i * dx; + double neighbor_right_x = current_x + dx, + neighbor_left_x = current_x - dx; + + if (neighbor_left_x < p.x && p.x <= current_x) + return Point{ (p.x - neighbor_left_x) / dx, 0 }; + + if (current_x < p.x && p.x < neighbor_right_x) + return Point{ (neighbor_right_x - p.x) / dx, 0 }; + + return Point{ 0, 0 }; + }); + + fem.set_AF(A, F); + fem.apply_boundaries(); + basis[2 * (2 * n_side_y + n_side_x + i)] = fem.solve(); + + // находим базисную функцию для, равную 1 в (2 * n_side_y + n_side_x + i) по y компоненте + fem.set_boundaries('W', zero); + fem.set_boundaries('N', zero); + fem.set_boundaries('E', zero); + fem.set_boundaries('S', + [&](const Point& p) + { + double current_x = b.x - i * dx; + double neighbor_right_x = current_x + dx, + neighbor_left_x = current_x - dx; + + if (neighbor_left_x < p.x && p.x <= current_x) + return Point{ 0, (p.x - neighbor_left_x) / dx }; + + if (current_x < p.x && p.x < neighbor_right_x) + return Point{ 0, (neighbor_right_x - p.x) / dx }; + + return Point{ 0, 0 }; + }); + + fem.set_AF(A, F); + fem.apply_boundaries(); + basis[2 * (2 * n_side_y + n_side_x + i) + 1] = fem.solve(); + + } + /*for (size_t i = 0; i < basis.size(); ++i){ + for (size_t j = 0; j < basis[0].size(); ++j) + std::cout << basis[i][j] << ' '; + std::cout << '\n'; + }*/ + Matrix W = matrix_form_basis(); + //W.print(); + //std::cout << "\n\n"; + K = W.T().dot(A).dot(W); +} + +Matrix FSEM::matrix_form_basis() { + Matrix W(2 * basis[0].size(), basis.size()); + for (size_t i = 0; i < basis[0].size(); ++i) + for (size_t j = 0; j < basis.size(); ++j) { + W[2 * i][j] = basis[j][i].x; + W[2 * i + 1][j] = basis[j][i].y; + } + return W; +} + +void FSEM::construct_f_bc2(const std::vector& pos, + const std::vector& g) { + Matrix W = matrix_form_basis(); + std::vector p_vec(W.size(), 0.0); + fem.calculate_bc2(pos, g, p_vec); + f = W.T().dot(p_vec); +} + +void FSEM::set_bc1(char side, const vec_function& g) { + if (side == 'W') + for (size_t i = 0; i <= n_side_y; ++i) + basis_coefficients[i] = g(nodes[i]); + + else if (side == 'N') + for (size_t i = n_side_y; i <= n_side_x + n_side_y; ++i) + basis_coefficients[i] = g(nodes[i]); + + else if (side == 'E') + for (size_t i = n_side_x + n_side_y; i <= n_side_x + 2 * n_side_y; ++i) + basis_coefficients[i] = g(nodes[i]); + + else { + for (size_t i = n_side_x + 2 * n_side_y; i < 2 * n_side_x + 2 * n_side_y; ++i) + basis_coefficients[i] = g(nodes[i]); + basis_coefficients[0] = g(nodes[0]); + } +} + +void FSEM::calculate_coef_Matrix_bc2(const int finish, const int i, + bool is_cur_Neumann, bool is_prev_Neumann, int& add_N, int& add_D, const int add_basis, Matrix& N, + Matrix& D) { + + // обработка первой точки на данной границе + if (is_cur_Neumann && is_prev_Neumann) { + N[2 * i][add_N] = basis[add_basis][i].x; + N[2 * i + 1][add_N] = basis[add_basis][i].y; + + N[2 * i][1 + add_N] = basis[1 + add_basis][i].x; + N[2 * i + 1][1 + add_N] = basis[1 + add_basis][i].y; + } + else { + D[2 * i][add_D] = basis[add_basis][i].x; + D[2 * i + 1][add_D] = basis[add_basis][i].y; + + D[2 * i][1 + add_D] = basis[1 + add_basis][i].x; + D[2 * i + 1][1 + add_D] = basis[1 + add_basis][i].y; + } + + if (is_cur_Neumann && !is_prev_Neumann) + add_N -= 2; + + for (int j = 2; j < finish; j++) { + if (is_cur_Neumann) { + N[2 * i][j + add_N] = basis[j + add_basis][i].x; + N[2 * i + 1][j + add_N] = basis[j + add_basis][i].y; + } + else { + D[2 * i][j + add_D] = basis[j + add_basis][i].x; + D[2 * i + 1][j + add_D] = basis[j + add_basis][i].y; + } + } + + if (is_cur_Neumann) { + add_N += finish; + if (!is_prev_Neumann) { + add_D += 2; + } + } + else + add_D += finish; +} + +void FSEM::save_bc1(std::vector& coefs_Dirichle, int& dir_id, const int i) { + coefs_Dirichle[dir_id] = basis_coefficients[i].x; + coefs_Dirichle[dir_id + 1] = basis_coefficients[i].y; + + dir_id += 2; +} + +void FSEM::set_bc2(const std::vector& pos, + const std::vector& g) { + + fem.construct_AF(E, nu, zero); + auto K_fem = fem.get_AF().first; // матрица жесткости + + size_t n_known_coefs = 0; + if (!pos[0]) n_known_coefs += n_side_y + 1; + if (!pos[2]) n_known_coefs += n_side_y + 1; + if (!pos[1]) n_known_coefs += n_side_x + 1; + if (!pos[3]) n_known_coefs += n_side_x + 1; + + if (!pos[0] && !pos[1]) n_known_coefs -= 1; + if (!pos[0] && !pos[3]) n_known_coefs -= 1; + if (!pos[2] && !pos[1]) n_known_coefs -= 1; + if (!pos[2] && !pos[3]) n_known_coefs -= 1; + + size_t n_unknown_coefs = 2 * (2 * n_side_x + 2 * n_side_y - n_known_coefs); + + // столбцы - значения суперэлементов, соответствующих + // неизвестным коэффициентам, в узлах мкэ сетки + Matrix N(K_fem.size(), n_unknown_coefs); + + // столбцы - значения суперэлементов, соответствующих + // известным коэффициентам, в узлах мкэ сетки + Matrix D(K_fem.size(), 2 * n_known_coefs); + + // интегралы от ГУ 2 рода * функции формы мкэ + std::vector p_vec(K_fem.size(), 0.0); + + int finish = 0, prev_pos = 0; + + // находим B и C + for (size_t i = 0; i < fem.psize(); i++) { + int add_N = 0, add_D = 0, add_basis = 0; + + for (size_t j = 0; j < 4; ++j) { + if (j == 0) + prev_pos = 3; + else + prev_pos = static_cast(j - 1); + + if (j == 0 || j == 2) + finish = static_cast(2 * n_side_y); + else + finish = static_cast(2 * n_side_x); + + calculate_coef_Matrix_bc2( + finish, + static_cast(i), + pos[j], + pos[prev_pos], + add_N, + add_D, + add_basis, + N, + D); + add_basis += finish; + } + } + + fem.calculate_bc2(pos, g, p_vec); + + Matrix N_Transposed = N.T(); + + Matrix A = N_Transposed.dot(K_fem).dot(N); + + // сохрвняем известные коэффициенты из ГУ Дирихле + std::vector coefs_Dirichle(2 * n_known_coefs); + int dir_id = 0; + + if (pos[0] && !pos[3]) + save_bc1(coefs_Dirichle, dir_id, 0); + + if (!pos[0]) { + for (int i = 0; i < n_side_y; ++i) + save_bc1(coefs_Dirichle, dir_id, i); + + if (pos[1]) + save_bc1(coefs_Dirichle, dir_id, static_cast(n_side_y)); + } + + if (!pos[1]) { + + for (int i = static_cast(n_side_y); + i < static_cast(n_side_y + n_side_x); ++i) + save_bc1(coefs_Dirichle, dir_id, i); + + if (pos[2]) + save_bc1(coefs_Dirichle, dir_id, static_cast(n_side_y + n_side_x)); + } + + if (!pos[2]) { + for (int i = static_cast(n_side_x + n_side_y); + i < static_cast(n_side_x + 2 * n_side_y); ++i) + save_bc1(coefs_Dirichle, dir_id, i); + + if (pos[3]) + save_bc1(coefs_Dirichle, dir_id, static_cast(n_side_x + 2 * n_side_y)); + } + + if (!pos[3]) + for (int i = static_cast(n_side_x + 2 * n_side_y); + i < static_cast(2 * (n_side_x + n_side_y)); ++i) + save_bc1(coefs_Dirichle, dir_id, i); + + std::vector f_fem; + if (pos[0] + pos[1] + pos[2] + pos[3] == 4) + f_fem = N_Transposed.dot(p_vec); + else + f_fem = N_Transposed.dot(p_vec) - N_Transposed.dot(K_fem).dot(D).dot(coefs_Dirichle); + + auto [L, U] = LU_decomposition(A); + std::vector ans = solveLU(L, U, f_fem); + + int ans_id = 0; + if (pos[0]) + for (size_t i = 0; i <= n_side_y; ++i) { + if (i == 0 && !pos[3]) + continue; + if (i == n_side_y && !pos[1]) + continue; + basis_coefficients[i] = { ans[ans_id], ans[ans_id + 1] }; + ans_id += 2; + } + + if (pos[1]) + for (size_t i = n_side_y + 1; i <= n_side_x + n_side_y; ++i) { + if (i == n_side_x + n_side_y && !pos[2]) + continue; + + basis_coefficients[i] = { ans[ans_id], ans[ans_id + 1] }; + ans_id += 2; + } + + if (pos[2]) + for (size_t i = n_side_x + n_side_y + 1; i <= n_side_x + 2 * n_side_y; ++i) { + if (i == n_side_x + 2 * n_side_y && !pos[3]) + continue; + basis_coefficients[i] = { ans[ans_id], ans[ans_id + 1] }; + ans_id += 2; + } + + if (pos[3]) { + for (size_t i = n_side_x + 2 * n_side_y + 1; i < 2 * n_side_x + 2 * n_side_y; ++i) { + basis_coefficients[i] = { ans[ans_id], ans[ans_id + 1] }; + ans_id += 2; + } + } + +} + +std::vector FSEM::find_answer() { + + std::vector res(fem.psize()); + for (size_t i = 0; i < res.size(); i++) + for (size_t j = 0; j < nodes.size(); j++) { + res[i].x += basis_coefficients[j].x * basis[2 * j][i].x + + basis_coefficients[j].y * basis[2 * j + 1][i].x; + res[i].y += basis_coefficients[j].x * basis[2 * j][i].y + + basis_coefficients[j].y * basis[2 * j + 1][i].y; + } + + return res; +} + +std::vector FSEM::get_side_nodes(char side) const { + std::vector side_nodes; + double value = 0; + bool hor = false; + + if (side == 'N') { + hor = true; + value = b.y; + } + else if (side == 'S') { + hor = true; + value = a.y; + } + else if (side == 'W') + value = a.x; + + else if (side == 'E') + value = b.x; + + for (size_t i = 0; i < nodes.size(); ++i) { + if (hor) { + side_nodes.reserve(n_side_x + 1); + if (fabs(nodes[i].y - value) < 1e-10) + side_nodes.push_back(i); + } + else { + side_nodes.reserve(n_side_y + 1); + if (fabs(nodes[i].x - value) < 1e-10) + side_nodes.push_back(i); + } + } + + std::sort(side_nodes.begin(), side_nodes.end(), [&](size_t lhs, size_t rhs) { + if (hor) + return nodes[lhs].x < nodes[rhs].x; + return nodes[lhs].y < nodes[rhs].y; + }); + + return side_nodes; +} + +std::vector FSEM::get_side_fem_nodes(char side) const { + std::vector side_nodes; + const size_t mx = coef_x * n_side_x + 1; + const size_t ny = coef_y * n_side_y + 1; + + if (side == 'S') { + side_nodes.reserve(n_side_x + 1); + for (size_t j = 0; j <= n_side_x; ++j) + side_nodes.push_back(j * coef_x); + } + else if (side == 'N') { + side_nodes.reserve(n_side_x + 1); + for (size_t j = 0; j <= n_side_x; ++j) + side_nodes.push_back(mx * (ny - 1) + j * coef_x); + } + else if (side == 'W') { + side_nodes.reserve(n_side_y + 1); + for (size_t i = 0; i <= n_side_y; ++i) + side_nodes.push_back(i * coef_y * mx); + } + else if (side == 'E') { + side_nodes.reserve(n_side_y + 1); + for (size_t i = 0; i <= n_side_y; ++i) + side_nodes.push_back(i * coef_y * mx + mx - 1); + } + + return side_nodes; +} + +std::vector> FSEM::get_known_dofs() const { + std::vector> known; + known.reserve(2 * basis_coefficients.size()); + + for (size_t i = 0; i < basis_coefficients.size(); ++i) { + if (!std::isnan(basis_coefficients[i].x)) + known.emplace_back(2 * i, basis_coefficients[i].x); + if (!std::isnan(basis_coefficients[i].y)) + known.emplace_back(2 * i + 1, basis_coefficients[i].y); + } + + return known; +} + +double mortar_shape_func(size_t i, const std::vector& s, double cur) { + if (i > 0 && cur >= s[i - 1] && cur <= s[i]) + return (cur - s[i - 1]) / (s[i] - s[i - 1]); + + if (i + 1 < s.size() && cur >= s[i] && cur <= s[i + 1]) + return (s[i + 1] - cur) / (s[i + 1] - s[i]); + return 0; +} + +namespace { + + constexpr double kContactEps = 1e-12; + constexpr double kAxisymmetricTwoPi = 6.28318530717958647692; + + size_t find_segment_index(const std::vector& x_nodes, double x_mid) { + if (x_nodes.size() < 2) + throw std::runtime_error("Contact boundary has less than two nodes."); + + for (size_t i = 0; i + 1 < x_nodes.size(); ++i) { + if (x_mid >= x_nodes[i] - kContactEps && x_mid <= x_nodes[i + 1] + kContactEps) + return i; + } + + throw std::runtime_error("Contact quadrature point is outside boundary segmentation."); + } + + double interpolate_trace_value( + const std::vector& basis_component, + const std::vector& fem_side_nodes, + const std::vector& side_x, + double x) { + + const size_t segment = find_segment_index(side_x, x); + const double x_left = side_x[segment]; + const double x_right = side_x[segment + 1]; + + const size_t left_fem = fem_side_nodes[segment]; + const size_t right_fem = fem_side_nodes[segment + 1]; + const double t = (x - x_left) / (x_right - x_left); + + return (1.0 - t) * basis_component[left_fem].y + + t * basis_component[right_fem].y; + } + + void assemble_body_mortar_matrix( + Matrix& M, + const std::vector>& basis, + const std::vector& side_nodes, + const std::vector& fem_side_nodes, + const std::vector& side_x, + const std::vector& mortar_elements, + const std::vector& lambda_nodes, + CoordinateSystem coordinate_system) { + + for (const auto& mortar_element : mortar_elements) { + const double r_left = mortar_element.chi_left; + const double r_right = mortar_element.chi_right; + const double len = r_right - r_left; + if (len <= kContactEps) + continue; + + const double r_mid = 0.5 * (r_left + r_right); + std::vector lambda_values(lambda_nodes.size()); + for (size_t l = 0; l < lambda_nodes.size(); ++l) + lambda_values[l] = mortar_shape_func(l, lambda_nodes, r_mid); + + for (size_t node : side_nodes) { + const double N_left_x = interpolate_trace_value( + basis[2 * node], fem_side_nodes, side_x, r_left); + const double N_right_x = interpolate_trace_value( + basis[2 * node], fem_side_nodes, side_x, r_right); + const double N_left_y = interpolate_trace_value( + basis[2 * node + 1], fem_side_nodes, side_x, r_left); + const double N_right_y = interpolate_trace_value( + basis[2 * node + 1], fem_side_nodes, side_x, r_right); + const double N_val_x = 0.5 * (N_left_x + N_right_x); + const double N_val_y = 0.5 * (N_left_y + N_right_y); + const double weight = (coordinate_system == CoordinateSystem::Axisymmetric) + ? kAxisymmetricTwoPi * r_mid * len + : len; + + for (size_t l = 0; l < lambda_nodes.size(); ++l) { + M[2 * node][l] += N_val_x * lambda_values[l] * weight; + M[2 * node + 1][l] += N_val_y * lambda_values[l] * weight; + } + } + } + } + +} + +std::vector solve_mortar_contact( + FSEM& bottom_body, + FSEM& top_body, + const std::vector& rhs_bottom, + const std::vector& rhs_top, + const ContactOptions& options) { + + Matrix A1 = bottom_body.get_K(); + Matrix A2 = top_body.get_K(); + + const auto& basis_bottom = bottom_body.get_basis(); + const auto& basis_top = top_body.get_basis(); + + std::vector side_bottom = bottom_body.get_side_nodes('N'); + std::vector side_top = top_body.get_side_nodes('S'); + std::vector fem_bottom = bottom_body.get_side_fem_nodes('N'); + std::vector fem_top = top_body.get_side_fem_nodes('S'); + + const size_t n1 = A1.size(); + const size_t n2 = A2.size(); + + std::vector bottom_x(side_bottom.size()); + std::vector top_x(side_top.size()); + for (size_t i = 0; i < side_bottom.size(); ++i) + bottom_x[i] = bottom_body[side_bottom[i]].x; + for (size_t i = 0; i < side_top.size(); ++i) + top_x[i] = top_body[side_top[i]].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()); + if (contact_left > contact_right + kContactEps) + throw std::runtime_error("Bodies do not overlap along the contact boundary."); + + 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); + } + + const std::vector& lambda_nodes = discretization.lambda_nodes; + const std::vector& mortar_elements = discretization.mortar_elements; + const size_t n_lambda = lambda_nodes.size(); + + const auto known_bottom = bottom_body.get_known_dofs(); + const auto known_top = top_body.get_known_dofs(); + + Matrix M1(n1, n_lambda); + Matrix M2(n2, n_lambda); + + assemble_body_mortar_matrix(M1, basis_bottom, side_bottom, fem_bottom, + bottom_x, mortar_elements, lambda_nodes, options.coordinate_system); + assemble_body_mortar_matrix(M2, basis_top, side_top, fem_top, + top_x, mortar_elements, lambda_nodes, options.coordinate_system); + + const size_t total = n1 + n2 + n_lambda; + + Matrix Sys(total); + std::vector rhs(total, 0); + + for (size_t i = 0; i < n1; ++i) { + rhs[i] = rhs_bottom[i]; + for (size_t j = 0; j < n1; ++j) + Sys[i][j] = A1[i][j]; + } + + for (size_t i = 0; i < n2; ++i) { + rhs[n1 + i] = rhs_top[i]; + for (size_t j = 0; j < n2; ++j) + Sys[n1 + i][n1 + j] = A2[i][j]; + } + + for (size_t i = 0; i < n1; ++i) + for (size_t j = 0; j < n_lambda; ++j) { + Sys[i][n1 + n2 + j] = M1[i][j]; + Sys[n1 + n2 + j][i] = M1[i][j]; + } + + for (size_t i = 0; i < n2; ++i) + for (size_t j = 0; j < n_lambda; ++j) { + Sys[n1 + i][n1 + n2 + j] = -M2[i][j]; + Sys[n1 + n2 + j][n1 + i] = -M2[i][j]; + } + + auto apply_known_dof = [&](size_t dof, double value) { + for (size_t i = 0; i < Sys[0].size(); ++i) + Sys[dof][i] = 0; + + Sys[dof][dof] = 1; + rhs[dof] = value; + }; + + for (const auto& [dof, value] : known_bottom) + apply_known_dof(dof, value); + + for (const auto& [dof, value] : known_top) + apply_known_dof(n1 + dof, value); + + return solveWithLU(Sys, rhs); +} + +std::vector solve_mortar_contact( + FSEM& bottom_body, + FSEM& top_body, + const std::vector& rhs_bottom, + const std::vector& rhs_top, + size_t lambda_node_count) { + + ContactOptions options; + options.coordinate_system = bottom_body.get_coordinate_system(); + options.method = ContactMethod::UniformUnionPartition; + options.lambda_node_count = lambda_node_count; + return solve_mortar_contact(bottom_body, top_body, rhs_bottom, rhs_top, options); +} + +std::vector FSEM::find_answer(const std::vector& coefs, size_t start) { + + std::vector res(fem.psize()); + for (size_t i = 0; i < res.size(); i++) + for (size_t j = 0; j < nodes.size(); j++) { + Point coef = { coefs[start + 2 * j] , + coefs[start + 2 * j + 1] }; + res[i].x += coef.x * basis[2 * j][i].x + + coef.y * basis[2 * j + 1][i].x; + res[i].y += coef.x * basis[2 * j][i].y + + coef.y * basis[2 * j + 1][i].y; + } + + return res; +} diff --git a/src/LinearAlgebra.cpp b/src/LinearAlgebra.cpp new file mode 100644 index 0000000..1b6c93b --- /dev/null +++ b/src/LinearAlgebra.cpp @@ -0,0 +1,328 @@ +#include "Elasticity.h" + +#include + +Matrix::Matrix(size_t n, double a) { + matrix = std::vector>(n, std::vector(n, a)); +} + +Matrix::Matrix(size_t n, size_t m, double a) { + matrix = std::vector>(n, std::vector(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 operator-(const std::vector& a, + const std::vector& b) { + if (a.size() != b.size()) + std::cout << "Wrong vector size for subtraction\n"; + + std::vector 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 Matrix::dot(const std::vector&v) const { + if (v.size() != size(1)) + throw std::runtime_error("Wrong matrix size while multiplication!"); + + std::vector 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 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 solveLU(const Matrix& L, const Matrix& U, + const std::vector& b) { + std::vector 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 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 solveGaussFullPivot(const Matrix& A, + const std::vector& b, + double eps) { + + const size_t n = A.size(); + Matrix M(A); + std::vector rhs = b; + + std::vector 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 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 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& 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 solveLU_with_perm(const Matrix& L, + const Matrix& U, + const std::vector& row_perm, + const std::vector& b) { + const size_t n = L.size(); + std::vector rhs(n); + for (size_t i = 0; i < n; ++i) rhs[i] = b[row_perm[i]]; + + std::vector 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 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 solveWithLU(const Matrix& A, + const std::vector& b, + double eps) { + const size_t n = A.size(); + + Matrix L(n, n), U(n, n); + std::vector 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); +} diff --git a/src/main.cpp b/src/main.cpp new file mode 100644 index 0000000..884f3a3 --- /dev/null +++ b/src/main.cpp @@ -0,0 +1,775 @@ +#include "Elasticity.h" +#include "TestCases.h" +#include "ContactSlaveNodes.h" +#include "ContactUniformLambdaPartition.h" +#include "ContactUniformUnionPartition.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace { + +constexpr double kTraceEps = 1e-12; + +struct TestRunConfig { + std::string name; + std::optional bottom_x; + std::optional bottom_y; + std::optional top_x; + std::optional top_y; + std::optional lambda_nodes; +}; + +struct InputConfig { + CoordinateSystem coordinate_system = CoordinateSystem::Axisymmetric; + ContactMethod contact_method = ContactMethod::UniformUnionPartition; + ContactSlaveBody slave_body = ContactSlaveBody::Bottom; + std::vector tests; +}; + +std::string trim(const std::string& text) { + std::string value = text; + if (value.size() >= 3 && + static_cast(value[0]) == 0xEF && + static_cast(value[1]) == 0xBB && + static_cast(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(std::tolower(static_cast(ch))); + return value; +} + +std::vector split_list(std::string value) { + for (char& ch : value) + if (ch == ',' || ch == ';') + ch = ' '; + + std::vector 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& 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& 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& 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 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 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& 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& x_nodes, + const std::vector& 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 contact_output_grid( + const std::vector& bottom_x, + const std::vector& 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 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& 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 recover_side_normal_stress( + CoordinateSystem coordinate_system, + const FSEM& body, + const std::vector& 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 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 recover_side_normal_displacement( + const FSEM& body, + const std::vector& field, + char side) { + + const auto side_nodes = body.get_side_fem_nodes(side); + std::vector 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& 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& 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& bottom_field, + const FSEM& top, + const std::vector& top_field, + double E, + double nu) { + + const std::vector bottom_x = side_axis_coordinates(bottom, 'N', true); + const std::vector top_x = side_axis_coordinates(top, 'S', true); + const std::vector grid = contact_output_grid(bottom_x, top_x); + const std::vector bottom_sigma = + recover_side_normal_stress(coordinate_system, bottom, bottom_field, 'N', E, nu); + const std::vector 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& bottom_field, + const FSEM& top, + const std::vector& top_field) { + + const std::vector bottom_x = side_axis_coordinates(bottom, 'N', true); + const std::vector top_x = side_axis_coordinates(top, 'S', true); + const std::vector grid = contact_output_grid(bottom_x, top_x); + const std::vector bottom_u = recover_side_normal_displacement(bottom, bottom_field, 'N'); + const std::vector 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 lambda_nodes_for_output( + const FSEM& bottom, + const FSEM& top, + const ContactOptions& options) { + + const std::vector bottom_x = side_axis_coordinates(bottom, 'N', false); + const std::vector 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& solution, + size_t start, + const std::vector& 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& 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& 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 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 bottom_field = bottom.find_answer(solution); + const std::vector 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; +}