利用python进行有限元分析（三）Bar_element-coding in Python

import numpy as np
from scipy.linalg import eigh

def bar(num_elems):
    restrained_dofs = [0,]

    # element mass and stiffness matrices for a bar
    m = np.array([[2, 1], [1, 2]]) / (6. * num_elems)
    k = np.array([[1, -1], [-1, 1]]) * float(num_elems)

    # construct global mass and stiffness matrices
    M = np.zeros(shape=(num_elems + 1, num_elems + 1))
    K = np.zeros(shape=(num_elems + 1, num_elems + 1))

    # assembly of elements
    for i in range(num_elems):
        M_temp = np.zeros(shape=(num_elems + 1, num_elems + 1))
        K_temp = np.zeros(shape=(num_elems + 1, num_elems + 1))
        M_temp[i:i+2,i:i+2] = m
        K_temp[i:i+2,i:i+2] = k
        M += M_temp
        K += K_temp

    # remove the fixed degrees of freedom
    for dof in restrained_dofs:
        for i in [0, 1]:
            M = np.delete(M, dof, axis=i)
            K = np.delete(K, dof, axis=i)

    # eigenvalue problem
    eigenvalues, eigenvectors = eigh(K, M)
    frequencies = np.sqrt(eigenvalues)
    return M, K, frequencies, eigenvectors

M, K, frequencies, eigenvectors = bar(2)
print(frequencies)

输出结果与先前计算吻合