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.
126 lines
4.5 KiB
126 lines
4.5 KiB
# -*- coding: utf-8 -*-
|
|
"""
|
|
Created on Mon Nov 23 20:57:00 2020
|
|
|
|
@author: Liangchao Zhu
|
|
"""
|
|
import numpy as np
|
|
def shapeFunction():
|
|
# https://www.mathworks.com/matlabcentral/fileexchange/13508-multi-dimensional-gauss-points-and-weights
|
|
# https://www.mathworks.com/matlabcentral/fileexchange/6862-gauss3d
|
|
gp = 0.577350269189626
|
|
gausspoint = np.zeros((8, 3), dtype=np.float32)
|
|
for i in range(2):
|
|
for j in range(2):
|
|
for k in range(2):
|
|
gausspoint[i * 4 + j * 2 + k, :] = [(-1) ** i * gp, (-1) ** j * gp, (-1) ** k * gp]
|
|
wt = np.ones((8, 1), dtype=np.float32)
|
|
|
|
r = gausspoint[:, 0]
|
|
s = gausspoint[:, 1]
|
|
t = gausspoint[:, 2]
|
|
# shape functions, [N1,N2,N3,N4;] for every gauss point
|
|
shape = np.zeros((8, 8), dtype=np.float32)
|
|
shape[:, 0] = (1 - r) * (1 - s) * (1 - t) / 8
|
|
shape[:, 1] = (1 + r) * (1 - s) * (1 - t) / 8
|
|
shape[:, 2] = (1 - r) * (1 + s) * (1 - t) / 8
|
|
shape[:, 3] = (1 + r) * (1 + s) * (1 - t) / 8
|
|
|
|
shape[:, 4] = (1 - r) * (1 - s) * (1 + t) / 8
|
|
shape[:, 5] = (1 + r) * (1 - s) * (1 + t) / 8
|
|
shape[:, 6] = (1 - r) * (1 + s) * (1 + t) / 8
|
|
shape[:, 7] = (1 + r) * (1 + s) * (1 + t) / 8
|
|
|
|
# derivatives, [dN1dxi,dN2dxi,dN3dxi,dN4dxi;] for every gauss point
|
|
dshapedr = np.zeros((8, 8), dtype=np.float32)
|
|
dshapedr[:, 0] = -(1 - s) * (1 - t) / 8
|
|
dshapedr[:, 1] = (1 - s) * (1 - t) / 8
|
|
dshapedr[:, 2] = -(1 + s) * (1 - t) / 8
|
|
dshapedr[:, 3] = (1 + s) * (1 - t) / 8
|
|
|
|
dshapedr[:, 4] = -(1 - s) * (1 + t) / 8
|
|
dshapedr[:, 5] = (1 - s) * (1 + t) / 8
|
|
dshapedr[:, 6] = -(1 + s) * (1 + t) / 8
|
|
dshapedr[:, 7] = (1 + s) * (1 + t) / 8
|
|
|
|
# [dN1deta,dN2deta,dN3deta,dN4deta;] for every gauss point
|
|
dshapeds = np.zeros((8, 8), dtype=np.float32)
|
|
dshapeds[:, 0] = -(1 - r) * (1 - t) / 8
|
|
dshapeds[:, 1] = -(1 + r) * (1 - t) / 8
|
|
dshapeds[:, 2] = (1 - r) * (1 - t) / 8
|
|
dshapeds[:, 3] = (1 + r) * (1 - t) / 8
|
|
|
|
dshapeds[:, 4] = -(1 - r) * (1 + t) / 8
|
|
dshapeds[:, 5] = -(1 + r) * (1 + t) / 8
|
|
dshapeds[:, 6] = (1 - r) * (1 + t) / 8
|
|
dshapeds[:, 7] = (1 + r) * (1 + t) / 8
|
|
|
|
dshapedt = np.zeros((8, 8), dtype=np.float32)
|
|
dshapedt[:, 0] = -(1 - r) * (1 - s) / 8
|
|
dshapedt[:, 1] = -(1 + r) * (1 - s) / 8
|
|
dshapedt[:, 2] = -(1 - r) * (1 + s) / 8
|
|
dshapedt[:, 3] = -(1 + r) * (1 + s) / 8
|
|
|
|
dshapedt[:, 4] = (1 - r) * (1 - s) / 8
|
|
dshapedt[:, 5] = (1 + r) * (1 - s) / 8
|
|
dshapedt[:, 6] = (1 - r) * (1 + s) / 8
|
|
dshapedt[:, 7] = (1 + r) * (1 + s) / 8
|
|
return wt, shape, dshapedr, dshapeds, dshapedt
|
|
|
|
def ISOelasticitytensor(E = 1, nu = 0.33):
|
|
D = np.zeros((6, 6), dtype=np.float32)
|
|
D[0, 0] = 1 - nu
|
|
D[1, 1] = 1 - nu
|
|
D[2, 2] = 1 - nu
|
|
|
|
D[0, 1] = nu
|
|
D[0, 2] = nu
|
|
D[1, 0] = nu
|
|
D[1, 2] = nu
|
|
D[2, 0] = nu
|
|
D[2, 1] = nu
|
|
|
|
D[3, 3] = (1 - 2 * nu) / 2
|
|
D[4, 4] = (1 - 2 * nu) / 2
|
|
D[5, 5] = (1 - 2 * nu) / 2
|
|
|
|
D = D * E / ((1 + nu) * (1 - 2 * nu))
|
|
return D
|
|
def LocalIntegratedMatrices(D,B0, detjacob, wt, eleNum):
|
|
# print(D)
|
|
# K = np.zeros((eleNum, 24, 24), dtype=np.float32)
|
|
# F = np.zeros((eleNum, 24, 6), dtype=np.float32)
|
|
for i in range(eleNum):
|
|
Ke = np.zeros((24, 24), dtype=np.float32)
|
|
Fe = np.zeros((24, 6), dtype=np.float32)
|
|
intB = np.zeros((6,24), dtype=np.float32)
|
|
for j in range(8):
|
|
B0i = B0[i, j * 6 * 24:(j + 1) * 6 * 24]
|
|
B0i = np.reshape(B0i, (6, 24))
|
|
detJ = detjacob[i, j]
|
|
B0iT = B0i.transpose()
|
|
Ke += B0iT @ D @ B0i * detJ * wt[j]
|
|
Fe += B0iT @ D * detJ * wt[j]
|
|
intB += B0i * detJ * wt[j]
|
|
# K[i, :, :] = Ke
|
|
return Ke,Fe,intB
|
|
|
|
def LocalKeFe(resolution,D0):
|
|
h=1.0/resolution
|
|
nodes = np.array([[0,0,0],
|
|
[1,0,0],
|
|
[0,1,0],
|
|
[1,1,0],
|
|
[0,0,1],
|
|
[1,0,1],
|
|
[0,1,1],
|
|
[1,1,1]])*h
|
|
eles = np.array([[0,1,2,3,4,5,6,7]])
|
|
eleNum=1
|
|
from elasticity3Dhex import Jacobian,dNdx_dNdy_dNdz_G,kinematicStiffnessLinear
|
|
wt, shape, dshapedr, dshapeds, dshapedt = shapeFunction()
|
|
detjacob, invjacob = Jacobian(nodes, eles, dshapedr, dshapeds, dshapedt, eleNum)
|
|
dshapedx, dshapedy, dshapedz = dNdx_dNdy_dNdz_G(invjacob, dshapedr, dshapeds, dshapedt, eleNum)
|
|
B0 = kinematicStiffnessLinear(dshapedx, dshapedy, dshapedz,eleNum)
|
|
Ke,Fe,intB = LocalIntegratedMatrices(D0,B0, detjacob, wt, eleNum)
|
|
return Ke,Fe,intB
|