面向对象有限元编程|数值计算类_python

面向对象有限元编程|数值计算类

生成矩阵

python主要依赖第三方库numpy,其中np.array和np.mat有区别，主要体现在：

生成矩阵所需格式不同

np.mat可以从字符串或列表中生成，而np.array只能从列表中生成。


import numpy as np
a = np.mat("1,2; 3,4")  #字符串生成2x2矩阵
b = np.mat([ [5,6], [7,8] ])   #列表生成2x2矩阵

c = np.array([ [2,6], [5,8] ])

生成的数组运算方式不同

np.array生成矩阵，用np.dot()表示矩阵乘法，星号(*)或np.multiply()表示点乘(对应元素相乘)。


import numpy as np
a = np.array([ [1, 1],  [1, 1] ])
b = np.array([ [2, 2],  [2, 2] ])
c = np.dot(a, b)     #矩阵乘法
d = np.multiply(a, b)  #对应元素相乘
e = a*b     #对应元素相乘

np.mat生成矩阵，星号(*)和np.dot()表示矩阵乘法，np.multiply()表示点乘(对应元素相乘)。


import numpy as np
a = np.mat([[1, 1], [1, 1]])
b = np.mat([[2, 2], [2, 2]])

c = a*b     #矩阵乘法
d =np.dot(a, b)   #矩阵乘法
e = np.multiply(a, b)  #对应元素相乘

矩阵索引

用a[0][0]访问矩阵a中第一行第一列元素，注意索引值从0开始。

解方程组

用np.linalg.solve(A, b)解方程组Ax=b，例如

##方程组Ax=b的解
import numpy as np
A =np.array( [ [6, -1.5, 1], [-1.5, 1.875, -1.5], [1, -1.5, 2]   ])
b = np.array( [ 30, -10, 0 ] )
x = np.linalg.solve(A, b)
print("Solution", x)

C++的数值计算类可以用第三方数值计算库，比如Eigen, BLAS等。

来源：数值分析与有限元编程

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

单元对象是构成整个结构对象的基本要素，如杆单元，梁单元，板单元，壳单元等等。虽然单元形状和特性各不相同，但基本特征和功能是相同的。比如都具有一定的几何形状，通过节点与其它单元连接，包含材料信息，在结构分析中各单元皆以单元刚度矩阵的形式组装成整体结构。各种单元的层次关系如图所示▲单元层次关系# 抽象单元类(基类)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;}来源：数值分析与有限元编程

面向对象有限元编程|数值计算类

生成矩阵

矩阵索引

解方程组

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