面向对象有限元编程|模块化_python

面向对象有限元编程|模块化

当程序变长时，可以方便地将其拆分为多个文件以便更容易地执行维护。此外，还可以将实现一个功能的程序段单独放在一个文件中，而不将它们的定义复制到每个程序中。可以通过将计划经常使用的函数的定义放入文件，Python将其称为模块。然后可以导入模块中的定义进入其他模块、脚本或与通常导入的交互模式命令。

模块是一个包含Python定义和语句的文件。文件名是模块名，在模块中，模块的名称（字符串形式）可用作全局变量的值__name__模块的__name__属性:在作为程序运行的时候，__name__属性的值是’main’；而作为模块导入的时候，这个值就是该模块的名字了。在阅读代码时，你会发现很多代码中都有 if name==’main’ 这么一行代码语句，也就不难理解这句代码的意思了。

将节点类，单元类，结构类，材料类等等单独放到一个文件里作为模块，解决具体问题时，再导入。

这里TrussModule文件就是节点类，单元类，结构类，材料类。test文件是一个具体的问题求解。调用模块的功能时，需要在前面加上模块名字。比如

import math
math.cos(1.0)

自定义模块也是如此：

import TrussModule  # 导入模块

mat1 = TrussModule.Material( 29.5e4 ) # 材料对象

# 节点对象
nd1 = TrussModule.TrussNode2D( 1, 0,    0,   0,  0,  0,       0 )
nd2 = TrussModule.TrussNode2D( 2, 400,   0,   1,  0,  20000,   0 )
nd3 = TrussModule.TrussNode2D( 3, 400,  300,   1,  1,  0,       -25000 )
nd4 = TrussModule.TrussNode2D( 3, 0,   300,   0,  0,  0,       0 )

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

面向对象有限元编程|整体结构类(下)

结构对象的重要任务是组装结构刚度矩阵和整体节点力向量。而结构刚度矩阵和整体节点力向量又依赖一个重要的参数--结构总自由度或者总自由度(gdof)。计算结构总自由度时，可以先将约束(边界)排除，这种方法叫先处理法。也可以不排除边界条件计算结构自由度，这种方法叫后处理法。▲先处理法结构标识，有约束的地方自由度为0如何计算结构总自由度？以先处理法为例，在创建节点对象时，定义了成员变量RX和RY，分别表示节点X和Y方向的约束，若有约束则其值为0，若没有约束则其值为1.如图所示初始时对结构总自由度变量gdof赋值为0，对所有的节点遍历，若RX或者RY的值为1，则gdof自加1，与此同时将gdof赋值给RX或者RY。若RX或者RY的值为0，则gdof的值不变。如图所示 #计算gdof def getGdof(self): gdof = 0 for i in range( len(self.listNode) ): if self.listNode[i].RX == 1: gdof = gdof + 1 self.listNode[i].RX = gdof if self.listNode[i].RY == 1: gdof = gdof + 1 self.listNode[i].RY = gdof return gdof 如何组装总刚度矩阵？以下思路同样适用于组装整体节点力向量。定义一个一维数组elToStr，赋值为每个单元的节点自由度编号，这和单元位移向量v=[v1,θ1,v2,θ2]是一一对应的。比如下图中单元1的elToStr=[0,0,0,1],单元2的elToStr=[0,1,0,2],单元3的elToStr=[0,2,0,3]，它是由单元向结构转化的枢纽。gdof=3，则结构刚度矩阵就是3行3列。开始组装结构刚度矩阵时，结构刚度矩阵赋值全为0。单元1只有一个节点自由度θ2没有被约束，单元1的elToStr只有第四个元素不为0，借鉴划行划列思路，将单元刚度矩阵中被约束的部分划掉，通过数组elToStr将余下的元素添加到结构刚度矩阵的对应位置1行1列此时结构刚度矩阵为单元2有2个节点自由度θ1,θ2没有被约束，按照上述思路此时结构刚度矩阵为单元3有2个节点自由度θ1,θ2没有被约束，按照上述思路此时结构刚度矩阵为 #elToStr数组 def elemToStruct(self, i): elToStr = [0, 0, 0, 0] elToStr[0] = self.listElement[i].node1.RX elToStr[1] = self.listElement[i].node1.RY elToStr[2] = self.listElement[i].node2.RX elToStr[3] = self.listElement[i].node2.RY return elToStr#组装结构刚度矩阵KK def packageStiffnessMatrix(self, elToStr, ek, KK): for m in range ( len(elToStr) ): for n in range ( len(elToStr) ): if (elToStr[m] != 0) and (elToStr[n] != 0) : KK[ elToStr[m]-1 ][ elToStr[n]-1 ] = \ KK[ elToStr[m]-1 ][ elToStr[n]-1 ] + ek [ m ][ n ] 一个完整的桁架结构对象类class FEModel: def __init__ (self, listNode, listElement): self.listNode = listNode self.listElement = listElement self.FF = None def getGdof(self): gdof = 0 for i in range( len(self.listNode) ): if self.listNode[i].RX == 1: gdof = gdof + 1 self.listNode[i].RX = gdof if self.listNode[i].RY == 1: gdof = gdof + 1 self.listNode[i].RY = gdof return gdof #组装节点力向量 def packageForce(self): for i in range( len(self.listNode) ): if self.listNode[i].RX != 0: self.FF[self.listNode[i].RX-1 ] = self.listNode[i].fx if self.listNode[i].RY != 0 : self.FF[self.listNode[i].RY-1 ] = self.listNode[i].fy def elemToStruct(self, i): elToStr = [0, 0, 0, 0] elToStr[0] = self.listElement[i].node1.RX elToStr[1] = self.listElement[i].node1.RY elToStr[2] = self.listElement[i].node2.RX elToStr[3] = self.listElement[i].node2.RY return elToStr def packageStiffnessMatrix(self, i,elToStr, ek, KK): for m in range ( len(elToStr) ): for n in range ( len(elToStr) ): if (elToStr[m] != 0) and (elToStr[n] != 0) : KK[ elToStr[m]-1 ][ elToStr[n]-1 ] = \ KK[ elToStr[m]-1 ][ elToStr[n]-1 ] + ek [ m ][ n ] 来源：数值分析与有限元编程