You can not select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
209 lines
6.9 KiB
209 lines
6.9 KiB
#
|
|
# 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()
|