c++编写五节点四面体有限元法求结点位移

五节点四面体有限元法是一种常用的有限元分析方法，可以用来求解结构的位移、应力等物理量。下面是一个简单的C++代码，用于实现五节点四面体有限元法求结点位移：

#include <iostream>
#include <cmath>

using namespace std;

const double PI = 3.141592653589793;
const int NODE_NUM = 5;
const int ELEM_NUM = 1;

double E = 2.0e11;     // 弹性模量
double nu = 0.3;       // 泊松比
double rho = 7800.0;   // 密度
double g = 9.81;       // 重力加速度

double X[NODE_NUM] = {0.0, 1.0, 0.0, 0.0, 0.0};
double Y[NODE_NUM] = {0.0, 0.0, 1.0, 0.0, 0.0};
double Z[NODE_NUM] = {0.0, 0.0, 0.0, 1.0, 0.0};
int connect[ELEM_NUM][NODE_NUM] = { {0, 1, 2, 3, 4} };

double K[NODE_NUM * 3][NODE_NUM * 3];   // 系数矩阵
double F[NODE_NUM * 3];                 // 右端向量
double U[NODE_NUM * 3];                 // 位移向量

double det4(double a[4][4])
{
    double b[3][3];
    for(int i=0; i<3; ++i)
        for(int j=0; j<3; ++j)
            b[i][j] = a[i+1][j+1];

    double cof1 = a[1][1]*a[2][2]*a[3][3] + a[1][2]*a[2][3]*a[3][1] + a[1][3]*a[2][1]*a[3][2];
    double cof2 = a[1][3]*a[2][2]*a[3][1] + a[1][1]*a[2][3]*a[3][2] + a[1][2]*a[2][1]*a[3][3];

    return cof1 - cof2;
}

void assemble()
{
    for(int i=0; i<NODE_NUM*3; ++i)
        for(int j=0; j<NODE_NUM*3; ++j)
            K[i][j] = 0.0;

    for(int e=0; e<ELEM_NUM; ++e)
    {
        int node[NODE_NUM];
        for(int i=0; i<NODE_NUM; ++i)
            node[i] = connect[e][i];

        double x[NODE_NUM], y[NODE_NUM], z[NODE_NUM];
        for(int i=0; i<NODE_NUM; ++i)
        {
            x[i] = X[node[i]];
            y[i] = Y[node[i]];
            z[i] = Z[node[i]];
        }

        double V = fabs(det4(  {{1.0, x[0], y[0], z[0]},
                               {1.0, x[1], y[1], z[1]},
                               {1.0, x[2], y[2], z[2]},
                               {1.0, x[3], y[3], z[3]}} ));

        // 计算单元刚度矩阵
        double Ke[NODE_NUM*3][NODE_NUM*3];
        for(int i=0; i<NODE_NUM*3; ++i)
            for(int j=0; j<NODE_NUM*3; ++j)
                Ke[i][j] = 0.0;

        double C = E / (1.0 - nu*nu);
        double D = C * nu;
        double E0 = C * (1.0 - nu) / (1.0 + nu) / (1.0 - 2.0*nu);
        double E1 = C / 2.0 / (1.0 + nu);

        double B[6][NODE_NUM*3];
        for(int i=0; i<6; ++i)
            for(int j=0; j<NODE_NUM*3; ++j)
                B[i][j] = 0.0;

        B[0][0] = B[3][3] = B[5][4] = 1.0 / V;
        B[1][1] = B[4][4] = B[5][3] = 1.0 / V;
        B[2][2] = B[4][3] = B[5][2] = 1.0 / V;
        B[3][1] = B[4][0] = B[5][1] = -1.0 / V;
        B[4][2] = B[5][0] = B[0][3] = -1.0 / V;
        B[5][4] = B[0][2] = B[1][0] = -1.0 / V;

        for(int i=0; i<NODE_NUM; ++i)
        {
            int ix = i * 3;
            Ke[ix][ix] += C * V;
            Ke[ix+1][ix+1] += C * V;
            Ke[ix+2][ix+2] += C * V;
        }

        for(int i=0; i<6; ++i)
        {
            double tmp[NODE_NUM*3];
            for(int j=0; j<NODE_NUM*3; ++j)
                tmp[j] = 0.0;

            for(int j=0; j<NODE_NUM; ++j)
            {
                int ix = j * 3;
                tmp[ix] = B[i][ix];
                tmp[ix+1] = B[i][ix+1];
                tmp[ix+2] = B[i][ix+2];
            }

            for(int j=0; j<NODE_NUM*3; ++j)
                for(int k=0; k<NODE_NUM*3; ++k)
                    Ke[j][k] += E0 * tmp[j] * tmp[k] * V + E1 * tmp[j] * tmp[k];
        }

        // 将单元刚度矩阵组装到全局刚度矩阵中
        for(int i=0; i<NODE_NUM; ++i)
        {
            int ix = node[i] * 3;
            for(int j=0; j<NODE_NUM; ++j)
            {
                int jx = node[j] * 3;
                for(int k=0; k<3; ++k)
                {
                    int ijk = ix + k;
                    int jkj = jx + k;
                    K[ijk][jkj] += Ke[ijk][jkj];
                }
            }
        }
    }
}

void solve()
{
    // 初始化右端向量
    for(int i=0; i<NODE_NUM*3; ++i)
        F[i] = 0.0;

    // 在右端向量中加入重力荷载
    for(int i=0; i<NODE_NUM; ++i)
    {
        int ix = i * 3;
        F[ix+1] -= rho * g;
    }

    // 固定边界条件
    for(int i=0; i<NODE_NUM; ++i)
    {
        int ix = i * 3;
        if(X[i] == 0.0 && Y[i] == 0.0 && Z[i] == 0.0)
        {
            for(int j=0; j<NODE_NUM*3; ++j)
            {
                K[ix][j] = 0.0;
                K[ix+1][j] = 0.0;
                K[ix+2][j] = 0.0;
            }
            K[ix][ix] = 1.0;
            K[ix+1][ix+1] = 1.0;
            K[ix+2][ix+2] = 1.0;
            F[ix] = 0.0;
            F[ix+1] = 0.0;
            F[ix+2] = 0.0;
        }
    }

    // 求解位移向量
    for(int i=0; i<NODE_NUM*3; ++i)
        U[i] = 0.0;

    for(int i=0; i<NODE_NUM*3; ++i)
    {
        if(K[i][i] == 0.0)
            continue;

        for(int j=i+1; j<NODE_NUM*3; ++j)
        {
            double factor = K[j][i] / K[i][i];
            for(int k=i+1; k<NODE_NUM*3; ++k)
                K[j][k] -= factor * K[i][k];
            F[j] -= factor * F[i];
        }
    }

    for(int i=NODE_NUM*3-1; i>=0; --i)
    {
        if(K[i][i] == 0.0)
            continue;

        U[i] = F[i] / K[i][i];

        for(int j=i-1; j>=0; --j)
            F[j] -= K[j][i] * U[i];
    }
}

void print_result()
{
    cout << "Node     X-Disp     Y-Disp     Z-Disp" << endl;
    for(int i=0; i<NODE_NUM; ++i)
    {
        int ix = i * 3;
        printf("%4d %10.4f %10.4f %10.4f\n", i, U[ix], U[ix+1], U[ix+2]);
    }
}

int main()
{
    assemble();
    solve();
    print_result();

    return 0;
}

在这个代码中，我们首先定义了一些常量和数组，包括弹性模量、泊松比、密度等物理量，以及节点坐标、单元连接关系等几何信息。然后，我们实现了一个 `assemble()` 函数，用于组装全局刚度矩阵。在这个函数中，我们首先遍历所有的单元，计算出每个单元的刚度矩阵，然后将它们组装到全局刚度矩阵中。在计算单元刚度矩阵的过程中，我们需要首先计算出单元体积，然后根据公式计算 B 矩阵和单元刚度矩阵。

接下来，我们实现了一个 `solve()` 函数，用于求解位移向量。在这个函数中，我们首先将右端向量初始化为零，并加入重力荷载。然后，我们对固定边界条件进行处理，将相关行和列的系数矩阵元素清零，并在右端向量中加入位移边界条件。最后，我们通过高斯消元法求解位移向量。

最后，我们实现了一个 `print_result()` 函数，用于输出计算结果。这个函数只是简单地打印每个节点的位移值。

当你运行这个程序时，它会输出每个节点的 X、Y、Z 三个方向的位移值。