有限元一阶四面体单元python编程（一）

#一阶四面体实体单元
#共4个节点，每个节点有且仅有 UX,UY,UZ三个自由度import numpy as np
from numpy import dot #矩阵乘法
from numpy.linalg import det # 求行列式
class Tetra3D(object): 
   def __init__(self,nodes, E, nu):
        self.nodes = nodes
        self.E = E
        self.nu = nu
        self.sixV = self.six_V()
        self.B = self.calc_B()
        self.D = self.calc_D()
        self.Ke = self.calc_Ke()
        
    def six_V(self): # 四面体体积的6倍
        m = np.ones((4,4),dtype = np.float)
        for i in range(4):
            m[i,1:] = self.nodes[i]
        return det(m)    def calc_B(self):
        R = (0,1,2,3) # 替代range(4)
        
        #Alpha = [ (-1)**i * det(np.delete(self.nodes,i,axis=0)) for i in R] #用不上        a = np.ones((4,3))
        for i in R:
            a[i,1:] = self.nodes[i,(1,2)]
        Belta = [ (-1)**(i+1) *det(np.delete(a,i,axis=0)) for i in R]        a = np.ones((4,3))
        for i in R:
            a[i,1:] = self.nodes[i,(0,2)]
        Gama = [ (-1)**i *det(np.delete(a,i,axis=0)) for i in R]
            
        a = np.ones((4,3))
        for i in R:
            a[i,1:] = self.nodes[i,(0,1)]
        Delta = [ (-1)**(i+1) *det(np.delete(a,i,axis=0)) for i in R]
        
        B = np.array(((Belta[0],    0,  0,  Belta[1],   0,  0,  Belta[2],   0,  0,  Belta[3],   0,  0),
                     (  0,      Gama[0],0,      0,  Gama[1],0,  0,      Gama[2],0,  0,  Gama[3],0),
                     ( 0,       0,  Delta[0],0,     0,      Delta[1],0,     0,Delta[2],0,0,Delta[3]),
                     (Gama[0],Belta[0],0,Gama[1],Belta[1],0,Gama[2],Belta[2],0,Gama[3],Belta[3],0,),
                     (0,Delta[0],Gama[0],0,Delta[1],Gama[1],0,Delta[2],Gama[2],0,Delta[3],Gama[3]),
                     (Delta[0],0,Belta[0],Delta[1],0,Belta[1],Delta[2],0,Belta[2],Delta[3],0,Belta[3])),
                     dtype = np.float64)
        B *= 1.0/self.sixV        return B
        
    def calc_D(self):
        nu = self.nu
        a = 1-nu
        b = (1-2*nu)/2.0
        D  = np.array(((a, nu, nu, 0, 0, 0),
                       (nu,a,  nu, 0, 0, 0),
                       (nu,nu, a,  0, 0, 0),
                       (0, 0,  0,  b, 0, 0),
                       (0, 0,  0,  0, b, 0),
                       (0, 0,  0,  0, 0, b)),
                       dtype = np.float64)
        D *= self.E/(1+nu)/(2*b)
        return D    def calc_Ke(self):
        Ke = self.sixV/6.0 * dot(self.B.T, dot(self.D,self.B))
        return Ke    def calc_Strain(self,disp):#disp大小12x1,单元位移矩阵
        return dot(self.B,disp)
    
    def calc_Stress(self,disp):
        Strain = self.calc_Strain(disp)
        return dot(self.D,Strain)

if __name__ =="__main__":
     
    # 节点坐标 x1,y1; x2,y2; x3,y3; x4,y4 ,须逆时针排    #import time
    #t0 = time.time()
    #t1 = time.time()    #print(f"耗时{t1-t0}")    nodes= np.array([[0.5,  0.5,0],[0.25, 0.25,0.4],[0,    0,0],[0.5,  0,0]])
    #nodes= np.array([[0,    0,0],[0.5,  0,0],[0.5,  0.5,0],[0.25, 0.25,0.4]])
    e = Tetra3D(nodes, E= 1.0E6, nu = 0.25)
    print(e.sixV)
    print(e.B[0,0])

from numpy import array, mat,zeros, double, integer,float64,sqrt,ravel
from numpy.linalg import solve #,det
from random import random
from tetra3d import Tetra3D
#from readFromFemap_Tetra import readNeuFile
import time
time1 = time.time()
#单位体系 N,m, kg

#Nodes info.: x,y,z. #Node Id must start from 1,the step is 1.
#读入节点坐标和单元信息
#NODE, ELEM = readNeuFile("3d.neu")
#NODE = array(NODE,dtype=double)
#ELEM = array(ELEM,dtype=integer)
dofs_per_node = 3 #x,y,z

NODE= array([[0,0,0],
           [0.025,0,0],
           [0,0.5,0],
           [0.025,0.5,0],
           [0.,0.,0.25],
           [0.025,0,0.25],
           [0,0.5,0.25],
           [0.025,0.5,0.25]],
          dtype=float64)
#Node ID 从0开始，步长1为1自增
node_qty = NODE.shape[0] #节点数量
#MAT ID,TYPE ID,4 Nodes' ID
ELEM = array([[1,1,3,5,0,1],
              [1,1,0,3,2,6],
              [1,1,5,4,6,0],
              [1,1,5,6,7,3],
              [1,1,0,5,3,6]],
         dtype=integer)
#单元ID 从0开始，步长1为1自增
elem_qty = ELEM.shape[0] #单元数量
#Boundary conditions

#节点自由度 dof CONSTRAINtrain
fixed_nodeID = [0,1,4,5]
CONSTRAIN = [(i,j,0.0) for i in fixed_nodeID for j in range(dofs_per_node)] #约束x,y,z
#外力 Force矩阵
F = mat(zeros((node_qty * dofs_per_node,1)), dtype=double)
F[dofs_per_node*2 + 1]= 3.125 # Fy on node 2,第1个自由度上（每节点3自由度时）
F[dofs_per_node*3 + 1]= 6.25 # Fy on node 3,第1个自由度上（每节点3自由度时）
F[dofs_per_node*7 + 1]= 3.125 # Fy on node 7,第1个自由度上（每节点3自由度时）
F[dofs_per_node*6 + 1]= 6.25 # Fy on node 6,第1个自由度上（每节点3自由度时）
print(F)
#材料属性
E= 210e6 # 弹性模量
nu = 0.30 #泊松比

time2 = time.time()
print(f"有限元模型输入完成！耗时{time2-time1}")
def assembleK(NODE,ELEM,CONSTRAIN):
    elems = [] #用于存储单元对象
    K= mat(zeros((node_qty*dofs_per_node ,node_qty *dofs_per_node)), dtype=float64)#初始化总刚度矩阵
    #遍历单元
    for ie in range(elem_qty):
        nodesID = ELEM[ie,2:6] #6 =2+4, 单元的4个节点的ID：
        i,j,m,n = nodesID
        nodes = array([NODE[i],NODE[j],NODE[m],NODE[n]])
        elem = Tetra3D(nodes,E, nu)#生成单元
        elems.append(elem)
        #总刚度矩阵组装（更新4x4个区域）
        M = 3
        for N1,I in enumerate(nodesID):
            for N2,J in enumerate(nodesID):
                K[dofs_per_node*I:dofs_per_node*I+M,dofs_per_node*J:dofs_per_node*J+M] += elem.Ke[M*N1:M*(N1+1), M*N2:M*(N2+1)]
        
    #将边界条件（力，位移约束）更新到刚度矩阵；更新外力矩阵
    BigNum = 1.0e150#大数法
    #cr = CONSTRAIN.shape[0]
    #a = (K-K.T)
    #print(K==K.T)
    
    #遍历约束节点
    for nodeID,dofID,nodeDisp in CONSTRAIN:
        jj = nodeID* dofs_per_node + dofID # 约束在总刚度矩阵的行号
        jj = int(jj)#CONSTRAIN 是浮点型数组，所有jj的结果是浮点型。做索引需是整数
        K[jj, jj] *= BigNum #总刚度矩阵对角线上的元素乘以大数
        F[jj] = K[jj, jj]*nodeDisp #载荷也用更新后的K[jj,jj]乘以指定位移
    return K, elems
       
# SOLVE
# 求解 位移矩阵Deformation Matrix DISP
#K,elems = assembleK(NODE,ELEM,CONSTRAIN)
K,elems= assembleK(NODE,ELEM,CONSTRAIN)
time3 = time.time()
print(f"总刚度矩阵组装完成！耗时{time3-time2}")

DISP= solve(K, F)#位移矩阵
#将位移在边界条件节点上的值用输入的约束值修正
for nodeID,dofID,nodeDisp in CONSTRAIN:
    jj = nodeID* dofs_per_node + dofID
    jj=int(jj)
    DISP[jj] = nodeDisp
time4 = time.time()
print(f"位移矩阵求解完成！耗时{time4-time3}")
print("位移矩阵：")
print(DISP)


#提高可读性：现在每个节点Ux,Uy，Uz 显示在同一行
Delta = DISP.reshape((-1,3))
X= NODE[:,0]
Y= NODE[:,1]
Z= NODE[:,2]
Delta_X = array(Delta[:,0])
Delta_Y = array(Delta[:,1])
Delta_Z = array(Delta[:,2])
#节点总位移 
DISP_total = sqrt(Delta_X**2 + Delta_Y**2,Delta_Z**2)
Stress_elems = zeros((elem_qty,6),dtype=float64)
#遍历单元,求各个单元上各个节点上的应力和应变
for ie in range(elem_qty):
    i,j,m,n = ELEM[ie,2:2+4] #单元的4个节点的ID：
    print(DISP[dofs_per_node*i+2,0])
    DISPe = array([[DISP[dofs_per_node*i,0]],
                [DISP[dofs_per_node*i+1,0]],
                [DISP[dofs_per_node*i+2,0]],
                [DISP[dofs_per_node*i+2,0]],
                [DISP[dofs_per_node*j+1,0]],
                [DISP[dofs_per_node*j+2,0]],
                [DISP[dofs_per_node*m,0]],
                [DISP[dofs_per_node*m+1,0]],
                [DISP[dofs_per_node*m+2,0]],
                [DISP[dofs_per_node*n,0]],
                [DISP[dofs_per_node*n+1,0]],
                [DISP[dofs_per_node*n+2,0]]])
    elem = elems[ie] #当前计算的单元

    #elem.calc_Strain(DISPe)
    Stress = elem.calc_Stress(DISPe)
    Stress_elems[ie] = ravel(Stress) #按C语言的顺序展平