用C++实现力密度法找形算法

admin1

1. 以下是使用C++实现力密度法找形算法的示例代码：
3. ```cpp
4. #include <iostream>
5. #include <vector>
6. #include <cmath>
8. using namespace std;
10. // 定义节点结构体
11. struct Node {
12.     double x, y; // 节点坐标
13.     double fx, fy; // 节点受力
14. };
16. // 计算两点之间的距离
17. double distance(Node a, Node b) {
18.     return sqrt((a.x - b.x) * (a.x - b.x) + (a.y - b.y) * (a.y - b.y));
19. }
21. // 计算节点受力
22. void calculateForces(vector<Node>& nodes, vector<vector<int>>& elements) {
23.     for (int i = 0; i < nodes.size(); i++) {
24.         nodes[i].fx = 0;
25.         nodes[i].fy = 0;
26.     }
28.     for (int i = 0; i < elements.size(); i++) {
29.         int n1 = elements[i][0];
30.         int n2 = elements[i][1];
31.         Node node1 = nodes[n1];
32.         Node node2 = nodes[n2];
33.         double L = distance(node1, node2);
34.         double P = 1000; // 假设施加在节点上的荷载为1000N/m
35.         double k = P / L; // 计算单元刚度
36.         double dx = node2.x - node1.x;
37.         double dy = node2.y - node1.y;
38.         double l = sqrt(dx * dx + dy * dy);
39.         double cosTheta = dx / l;
40.         double sinTheta = dy / l;
41.         double fx = k * (l - L) * cosTheta;
42.         double fy = k * (l - L) * sinTheta;
43.         nodes[n1].fx += fx;
44.         nodes[n1].fy += fy;
45.         nodes[n2].fx -= fx;
46.         nodes[n2].fy -= fy;
47.     }
48. }
50. // 求解线性方程组，获得自由节点坐标
51. void solveLinearSystem(vector<Node>& nodes, vector<vector<int>>& elements) {
52.     vector<vector<double>> A(nodes.size(), vector<double>(nodes.size()));
53.     vector<double> B(nodes.size());
55.     for (int i = 0; i < nodes.size(); i++) {
56.         A[i][i] = 1;
57.         B[i] = nodes[i].fx;
58.     }
60.     for (int i = 0; i < elements.size(); i++) {
61.         int n1 = elements[i][0];
62.         int n2 = elements[i][1];
63.         Node node1 = nodes[n1];
64.         Node node2 = nodes[n2];
65.         double L = distance(node1, node2);
66.         double P = 1000; // 假设施加在节点上的荷载为1000N/m
67.         double k = P / L; // 计算单元刚度
68.         double dx = node2.x - node1.x;
69.         double dy = node2.y - node1.y;
70.         double l = sqrt(dx * dx + dy * dy);
71.         double cosTheta = dx / l;
72.         double sinTheta = dy / l;
73.         double fx = k * (l - L) * cosTheta;
74.         double fy = k * (l - L) * sinTheta;
75.         A[n1][n1] += cosTheta * cosTheta;
76.         A[n1][n2] += cosTheta * sinTheta;
77.         A[n2][n1] += sinTheta * cosTheta;
78.         A[n2][n2] += sinTheta * sinTheta;
79.         B[n1] += fx * cosTheta;
80.         B[n2] += fx * sinTheta;
81.     }
83.     // 高斯消元法求解线性方程组
84.     for (int i = 0; i < nodes.size(); i++) {
85.         for (int j = i + 1; j < nodes.size(); j++) {
86.             double ratio = A[j][i] / A[i][i];
87.             B[j] -= ratio * B[i];
88.             for (int k = i; k < nodes.size(); k++) {
89.                 A[j][k] -= ratio * A[i][k];
90.             }
91.         }
92.     }
94.     for (int i = nodes.size() - 1; i >= 0; i--) {
95.         for (int j = i + 1; j < nodes.size(); j++) {
96.             B[i] -= A[i][j] * B[j];
97.         }
98.         B[i] /= A[i][i];
99.     }
101.     for (int i = 0; i < nodes.size(); i++) {
102.         nodes[i].x += B[i];
103.     }
104. }
106. int main() {
107.     // 定义节点和单元信息
108.     vector<Node> nodes = {{0, 0}, {5, 0}, {10, 0}, {15, 0}};
109.     vector<vector<int>> elements = {{0, 1}, {1, 2}, {2, 3}};
111.     // 计算节点受力
112.     calculateForces(nodes, elements);
114.     // 求解线性方程组，获得自由节点坐标
115.     solveLinearSystem(nodes, elements);
117.     // 输出结果
118.     for (int i = 0; i < nodes.size(); i++) {
119.         cout << "Node " << i << ": (" << nodes[i].x << ", " << nodes[i].y << ")" << endl;
120.     }
122.     return 0;
123. }
124. ```

复制代码