面向对象有限元编程|单元类

单元对象是构成整个结构对象的基本要素，如杆单元，梁单元，板单元，壳单元等等。虽然单元形状和特性各不相同，但基本特征和功能是相同的。比如都具有一定的几何形状，通过节点与其它单元连接，包含材料信息，在结构分析中各单元皆以单元刚度矩阵的形式组装成整体结构。各种单元的层次关系如图所示

▲单元层次关系

# 抽象单元类(基类)
class AbstractElement:
    def __init__ (self, id):
        self.id = id

对于一个桁架单元，其特征有编号、材料类型、横截面面积、两个节点的信息、单元长度、单元局部坐标与总体坐标的夹角、单元局部坐标与总体坐标的转换矩阵、单元局部坐标系下的刚度矩阵、单元整体坐标系下的刚度矩阵、应力矩阵、应变矩阵等等。

# 平面桁架单元类
class TrussElement2D(AbstractElement):
    def __init__ (self,id, A, mat, node1, node2):
        super().__init__(id)
        self.A = A
        self.mat = mat
        self.node1 = node1
        self.node2 = node2

    def elemDisp(self):
        return np.array( [ [0.], [0.], [0.], [0.] ] )
                          


    def elemLength(self):
        dx = self.node2.coord_X - self.node1.coord_X
        dy = self.node2.coord_Y - self.node1.coord_Y
        return  math.sqrt(dx*dx + dy*dy)

    def COS(self):
        dx = self.node2.coord_X - self.node1.coord_X
        return dx / self.elemLength()

    def SIN(self):
        dy = self.node2.coord_Y - self.node1.coord_Y
        return  dy / self.elemLength()

    def Trans(self):
        S = self.SIN()
        C = self.COS()
        return np.array([ [C, S, 0, 0],
                          [0, 0, C, S] ])

    def elemStiffnessMatrix(self):
        EA = self.mat.E * self.A
        S = self.SIN()
        C = self.COS()
        ek = np.array([ [C*C,  C*S, -C*C, -C*S ], 
                       [C*S,  S*S, -C*S, -S*S ], 
                       [-C*C, -C*S, C*C,  C*S], 
                       [-C*S, -S*S, C*S,  S*S]   ] )

        return EA/self.elemLength() * ek

    def strainMatrix(self):
        L  = self.elemLength()
        tmp1 = np.array( [-1.0, 1.0] )
        return 1/L * np.array( [-1.0, 1.0] )
        
    def normalForce(self, T, elemDisp):
        EA = self.mat.E * self.A
        b = self.strainMatrix() 
        tmp1 = EA * b
        
        tmp2 = np.dot(tmp1, T)
        return np.dot(tmp2, elemDisp)

采用从基类按层次继承来建立单元类的方法，除了重用父类代码之外，另一个很重要的优点就是可以运用多态。通过指向单元基类对象的指针，调用各派生单元类的成员函数。来看一个例子：从一个抽象单元类派生两个不同的单元类。现在要获取这两个单元的编号。

#include <iostream>
#include <vector>
using namespace std;

class AbstractElement
{
public:

    virtual ~AbstractElement() = default;
    virtual size_t getID() const = 0;   // 纯虚函数 
};


class Truss : public AbstractElement
{
protected:
    size_t id{ 0 };

public:
    Truss() = default;
    Truss(size_t id_) : id{ id_ } {}
    size_t getID() const override { return id; }


};


class Beam : public AbstractElement
{
protected:
    size_t id{ 0 };

public:
    Beam() = default;
    Beam(size_t id_) : id{ id_ } {}
    size_t getID() const override { return id; }
};

int main()
{
    Truss t1{ 1 };
    Beam  b1{ 2 };

    std::vector<AbstractElement*> elems{ &t1, &b1 };

    for (auto* el : elems)
    {
        std::cout << el->getID() << std::endl;  //这里体现了多态
    }

 
    return 0;

}