以二阶三节点单元为例，我们用形函数插值几何坐标：

对  求导：

因此微分关系是：

其中  是 Jacobian(1D为标量)

弱形式里的刚度矩阵需要 ，但形函数天然写在 ，所以要用链式法则把导数换回物理坐标：

3. 任意积分的统一计算

有了等参映射以后，任何  上的积分都能写成：

再用高斯积分得到计算公式：

4.选几个高斯点？

 点Gauss–Legendre在上对次数的多项式积分是精确的

实际 FEM 里高斯点个数的选择，取决于你要积分的被积函数最高有多复杂。

经验上可以这样取：

• • 刚度矩阵：若 ,  在单元内近似常数，线性单元用 1 点即可。
• • 体力项：如果  至少是线性的，线性单元通常需要 2 点才能算准。
• • 高阶单元一般取  较为合适。

2. 1D 有限元 Python 代码

有限元代码的核心任务是：

计算单元矩阵与向量 → 组装整体方程 → 施加边界条件 → 求解。

本 demo 以之前的 1D 杆问题为例，实现：

• • 自定义等分网格
• • 自定义单元阶次 
• • 计算出未知节点位移 
• • 生成位移场 及应力场 

1. 问题描述（problem.py）

仍取前面推文中的案例，具体参数：

• • 杆长 ， 为常数
• • 体力 （或 )
• • 左端 ：自然边界 
• • 右端 ：本质边界 

E = 100000.0
A = 1.0

def q(x):
    return 10.0 * A * x

natural_bcs = {0.0 : 10.0}
essential_bcs = {2.0 : 0.0}

2. 网格（mesh.py）

等分网格要输出两样东西：

• • coords: 全局节点坐标
• • conn: 每个单元包含哪些节点（连通性）

对  阶单元，单元间共享端点，所以总节点数：

import numpy as np
def make_mesh_1d(x0, xl, nel, p):
    """
    x0: float,起点坐标
    x1: float,终点坐标
    nel: int,单元数
    p: int,单元阶数
    """

    nnp = nel * p + 1
    coords = np.linspace(x0, xl, nnp)
    conn = np.zeros((nel, p+1), dtype=int)

    for e in range(nel):
        conn[e, :] = np.arange(e * p, e * p + (p + 1))

    return coords, conn

3. 单元形函数(shape1d.py)

等参坐标 ，等距节点：

Lagrange 形函数：

导数：

def lagrange_shape_1d(p, xi):
    """
    输入：
    p: int,单元阶数
    xi: 等参元坐标，在[-1,1]范围上
    ----------
    输出：
    N: (p+1, ) xi为输入时的形函数值
    dN_dxi: (p+1, ) xi为输入时形函数一阶导数值
    """

    N = np.ones(p+1, dtype=float)
    dN_dxi = np.zeros(p+1, dtype=float)
    xii = np.zeros(p+1, dtype=float)

    ### 对应公式Eq.9
    for i in range(p+1):
        xii[i] = -1 + 2 * i / p
    ### 对应公式Eq.10
    for i in range(p+1):
        for j in range(p+1):
            if i != j:
                N[i] = N[i] * (xi - xii[j]) / (xii[i] - xii[j])
    ### 对应公式Eq.11
    for i in range(p+1):
        term = 0
        for k in range(p+1):
            if k != i:
                prod = 1
                for j in range(p+1):
                    if (j != k) and (j != i):
                        prod = prod * (xi - xii[j]) / (xii[i] - xii[j])
                term += (1/(xii[i] - xii[k])) * prod
        dN_dxi[i] = term

    return N, dN_dxi

4. 数值积分（quad.py）

最小实现：1～3 高斯点，后续可补充。

def gauss_legendre(n_gauss):
    """
    输入：
    n_gauss: int,高斯积分点
    ----------
    输出：
    xi_g: (n_gauss, ) 对应的高斯点坐标
    w_g: (n_gauss, ) 对应的权重
    """
    xi_g = np.zeros(n_gauss)
    w_g = np.zeros(n_gauss)

    if n_gauss == 1:
        w_g[0] = 2.0
        xi_g[0] = 0.0
    elif n_gauss == 2:
        xi_g[:] = [-0.5773502691896257, 0.5773502691896257]
        w_g[:]  = [1.0, 1.0]
    elif n_gauss == 3:
        xi_g[:] = [-0.7745966692414834, 0.0, 0.7745966692414834]
        w_g[:]  = [0.5555555555555556, 0.8888888888888888, 0.5555555555555556]
    else:
        raise ValueError("n must be 1, 2, or 3 in this demo.")

    return xi_g, w_g

5. 单元矩阵与向量（element1d.py）

单元积分统一写成参考单元 上的求和：

每个高斯点的计算顺序很固定：
Step 1 形函数与导数：

Step 2 物理坐标：

Step 3 Jacobian：

Step 4  矩阵：

Step 5 累加 ：

Step 6 累加 ：

import numpy as np
from quad import *
from shape1d import *
from problem import *

def element_matrices(coords_e, p, n_gauss):
    """
    输入：
    coords_e: 单元物理坐标
    p: 单元阶次
    n_gauss: 高斯点个数
    ----------
    输出：
    Ke: (p+1, p+1) 单元刚度矩阵
    fe: (p+1, ) 单元外荷载向量
    """
    nen = p+1           # 单元节点数
    Ke = np.zeros((nen, nen))
    fe = np.zeros(nen)

    xi_g, w_g = gauss_legendre(n_gauss)

    for i in range(n_gauss):
        xi = xi_g[i]
        w = w_g[i]
        # Step 1: N, dN/dxi
        N, dN_dxi = lagrange_shape_1d(p, xi)
        # Step 2: x(xi)
        x = N @ coords_e
        # Step 3: J(xi)
        J = dN_dxi @ coords_e
        # Step 4: B = dN/dx
        dN_dx = dN_dxi / J
        B = dN_dx
        # Step 5: Ke
        Ke += (E*A) * np.outer(B, B) * (w * J)
        # Step 6: fe (body force)
        fe += N * q(x) * (w * J)

    return Ke, fe

6. 组装与求解（solver1d.py）

组装整体矩阵

前面刚好有一个变量 conn 定义是的各单元内节点的编号。