# # fem3d_tet_linear.py # # Created by Wei Chen on 8/12/21 # import numpy as np from scipy.sparse import coo_matrix from scipy.sparse.linalg import spsolve import math ''' fem3d_tet_linear: linear fem of 3d tetrahedron model @input: nvet: number of vertices nele: number of tetrahedron input/cor.txt: corrdinates of vertices input/elem.txt: indices of tetrahedron elements @output: output/u.txt: displacement of vertices output/cor1.txt: result corrdinates of vertices ''' def fem3d_tet_linear(): print('fem simulation start') eeps = 1e-8 # user-defined material properities E = 73.0 # Young's modulus nu = 0.17 # Poisson's ratio alpha = 5.5e-7 Tref = 293.15 # constitutive matrix D = E / (1. + nu) / (1. - 2. * nu) * np.array( [[1-nu, nu, nu, 0., 0., 0.], [nu, 1-nu, nu, 0., 0., 0.], [nu, nu, 1-nu, 0., 0., 0.], [0., 0., 0., (1.0-2.0*nu)/2.0, 0., 0.], [0., 0., 0., 0., (1.0-2.0*nu)/2.0, 0.], [0., 0., 0., 0., 0., (1.0-2.0*nu)/2.0]], dtype=np.float64) # read coordinates of vertices, tetrahedron information nvet, nele, cor, elem = readTet('input/cube.txt') T = readT(nvet, 'input/cube_T.txt') T = alpha * (T - Tref) # prepare finite element analysis ndof = nvet * 3 U = np.zeros((ndof), dtype=np.float64) F = np.zeros((ndof), dtype=np.float64) # user-defined fixed dofs fixednids = np.asarray(cor[:, 2]==0.0).nonzero()[0] fixeddofs = np.concatenate((fixednids*3, fixednids*3+1, fixednids*3+2), axis=0) freedofs = np.setdiff1d(np.arange(ndof), fixeddofs) print('compute DBC') # assembly the stiffness matrix edofMat = np.kron(3*elem, np.ones((1, 3))) + np.kron(np.tile(np.arange(3), 4), np.ones((nele, 1))) edofMat = edofMat.astype(np.int32) iK = np.kron(edofMat, np.ones((12, 1))).flatten() jK = np.kron(edofMat, np.ones((1, 12))).flatten() iK = iK.astype(np.int32) jK = jK.astype(np.int32) sK = np.empty((nele, 12*12), dtype=np.float64) for eleI in range(nele): elenids = np.array([elem[eleI, 0], elem[eleI, 1], elem[eleI, 2], elem[eleI, 3]], dtype=np.int32) eledofs = np.kron(3*elenids, np.ones((1, 3))) + np.tile(np.arange(3), 4) eledofs = eledofs.astype(np.int32) Xe = cor[elenids, :] Te = T[elenids] KE, Fe = initKE(D, Xe, Te) sK[eleI, :] = KE.flatten('F') F[eledofs] = F[eledofs] + Fe sK = sK.flatten() K = coo_matrix((sK, (iK, jK)), shape=(ndof, ndof)).tocsc() # FE-analysis print('solve') K = K[freedofs, :][:, freedofs] U[freedofs] = spsolve(K, F[freedofs]) # write result print('write displacement U.txt and new corrdinates cor1.txt') with open('output/U.txt', 'w') as f: for i in range(ndof): print(U[i], file=f) cor1 = cor + U.reshape(nvet, 3) with open('output/cor1.txt', 'w') as f: for vI in range(nvet): print(cor1[vI, 0], cor1[vI, 1], cor1[vI, 2], file=f) def readTet(filePath): with open(filePath, 'r') as f: line = f.readline() line_split = line.split() nvet = int(line_split[0]) line = f.readline() line_split = line.split() nele = int(line_split[0]) cor = np.empty((nvet, 3), dtype=np.float64) for vI in range(nvet): line = f.readline() line_split = line.split() cor[vI, 0] = float(line_split[0]) cor[vI, 1] = float(line_split[1]) cor[vI, 2] = float(line_split[2]) elem = np.empty((nele, 4), dtype=np.int32) for eleI in range(nele): line = f.readline() line_split = line.split() elem[eleI, 0] = int(line_split[0]) elem[eleI, 1] = int(line_split[1]) elem[eleI, 2] = int(line_split[2]) elem[eleI, 3] = int(line_split[3]) elem = elem - 1 return nvet, nele, cor, elem def readT(nvet, filePath): T = np.empty((nvet), dtype=np.float64) with open(filePath, 'r') as f: for vI in range(nvet): line = f.readline() line_split = line.split() T[vI] = float(line_split[0]) return T # generate element stiffness matrix # @input # D: constitutive matrix # X: 4*3 matrix, which is corridinates of element vertices # T: vector of 4 size # @output # KE: element stiffness matrix def initKE(D, X, T): H = np.array([[1., X[0, 0], X[0, 1], X[0, 2]], [1., X[1, 0], X[1, 1], X[1, 2]], [1., X[2, 0], X[2, 1], X[2, 2]], [1., X[3, 0], X[3, 1], X[3, 2]]]) V6 = abs(np.linalg.det(H)) a1 = np.linalg.det(H[[1, 2, 3], :][:, [1, 2, 3]]) a2 = -np.linalg.det(H[[0, 2, 3], :][:, [1, 2, 3]]) a3 = np.linalg.det(H[[0, 1, 3], :][:, [1, 2, 3]]) a4 = -np.linalg.det(H[[0, 1, 2], :][:, [1, 2, 3]]) b1 = -np.linalg.det(H[[1, 2, 3], :][:, [0, 2, 3]]) b2 = np.linalg.det(H[[0, 2, 3], :][:, [0, 2, 3]]) b3 = -np.linalg.det(H[[0, 1, 3], :][:, [0, 2, 3]]) b4 = np.linalg.det(H[[0, 1, 2], :][:, [0, 2, 3]]) c1 = np.linalg.det(H[[1, 2, 3], :][:, [0, 1, 3]]) c2 = -np.linalg.det(H[[0, 2, 3], :][:, [0, 1, 3]]) c3 = np.linalg.det(H[[0, 1, 3], :][:, [0, 1, 3]]) c4 = -np.linalg.det(H[[0, 1, 2], :][:, [0, 1, 3]]) d1 = -np.linalg.det(H[[1, 2, 3], :][:, [0, 1, 2]]) d2 = np.linalg.det(H[[0, 2, 3], :][:, [0, 1, 2]]) d3 = -np.linalg.det(H[[0, 1, 3], :][:, [0, 1, 2]]) d4 = np.linalg.det(H[[0, 1, 2], :][:, [0, 1, 2]]) B = np.empty((6, 12), dtype=np.float64) B[:, 0:3] = np.array([[b1, 0., 0.], [0., c1, 0.], [0., 0., d1], [c1, b1, 0.], [0., d1, c1], [d1, 0., b1]]) B[:, 3:6] = np.array([[b2, 0., 0.], [0., c2, 0.], [0., 0., d2], [c2, b2, 0.], [0., d2, c2], [d2, 0., b2]]) B[:, 6:9] = np.array([[b3, 0., 0.], [0., c3, 0.], [0., 0., d3], [c3, b3, 0.], [0., d3, c3], [d3, 0., b3]]) B[:, 9:12] = np.array([[b4, 0., 0.], [0., c4, 0.], [0., 0., d4], [c4, b4, 0.], [0., d4, c4], [d4, 0., b4]]) B = B / V6 KE = V6 / 6.0 * (B.T @ D @ B) strain = np.mean(T) * np.array([1., 1., 1., 0., 0., 0.]) Fe = V6 / 6.0 * (B.T @ D @ strain) return KE, Fe if __name__ == '__main__': fem3d_tet_linear()