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.
115 lines
3.2 KiB
115 lines
3.2 KiB
import numpy as np
|
|
from getDH import getDH
|
|
|
|
def getgp(s,d,n):
|
|
a=s
|
|
h=(d-s)/n
|
|
b=a+h
|
|
gp=np.zeros(n*2+1)
|
|
for i in range(1,n+1):
|
|
gp[i*2-1]=(a+b)/2+(b-a)/2/(-np.sqrt(3))
|
|
gp[i*2]=(a+b)/2+(b-a)/2/(np.sqrt(3))
|
|
a=b
|
|
b=b+h
|
|
return gp
|
|
|
|
|
|
def compute_dN(r,s,t):
|
|
dN=np.zeros((3+1,8+1))
|
|
#print(f's:{s}, shape: {type(s)}')
|
|
#dN[1,1]=s
|
|
dN[1,1]=-((s + 1)*(t - 1))/8
|
|
dN[1,2]=((s + 1)*(t - 1))/8
|
|
dN[1,3]=-((s - 1)*(t - 1))/8
|
|
dN[1,4]=((s - 1)*(t - 1))/8
|
|
dN[1,5]=((s + 1)*(t + 1))/8
|
|
dN[1,6]=-((s + 1)*(t + 1))/8
|
|
dN[1,7]=((s - 1)*(t + 1))/8
|
|
dN[1,8]=-((s - 1)*(t + 1))/8
|
|
dN[2,1]=-((r + 1)*(t - 1))/8
|
|
dN[2,2]=((r - 1)*(t - 1))/8
|
|
dN[2,3]=-((r - 1)*(t - 1))/8
|
|
dN[2,4]=((r + 1)*(t - 1))/8
|
|
dN[2,5]=((r + 1)*(t + 1))/8
|
|
dN[2,6]=-((r - 1)*(t + 1))/8
|
|
dN[2,7]=((r - 1)*(t + 1))/8
|
|
dN[2,8]=-((r + 1)*(t + 1))/8
|
|
dN[3,1]=-((r + 1)*(s + 1))/8
|
|
dN[3,2]=((r - 1)*(s + 1))/8
|
|
dN[3,3]=-((r - 1)*(s - 1))/8
|
|
dN[3,4]=((r + 1)*(s - 1))/8
|
|
dN[3,5]=((r + 1)*(s + 1))/8
|
|
dN[3,6]=-((r - 1)*(s + 1))/8
|
|
dN[3,7]=((r - 1)*(s - 1))/8
|
|
dN[3,8]=-((r + 1)*(s - 1))/8
|
|
return dN[1:,1:]
|
|
|
|
|
|
def LocalSMRT(X, Y, Z, C):
|
|
xs = -1
|
|
xd = 1
|
|
numgp_x = 2
|
|
gps_x = getgp(xs, xd, numgp_x)
|
|
K = np.zeros((24, 24))
|
|
for ix in range(1,numgp_x * 2+1):
|
|
x = gps_x[ix]
|
|
for iy in range(1,numgp_x * 2+1):
|
|
y = gps_x[iy]
|
|
for iz in range(1,numgp_x * 2+1):
|
|
z = gps_x[iz]
|
|
dN = compute_dN(x, y, z)
|
|
|
|
#dN=np.random.rand(3,8)
|
|
|
|
J = np.matmul(dN , np.asarray([X,Y,Z]).transpose())
|
|
#print (dN)
|
|
#print (np.asarray([X,Y,Z]).transpose())
|
|
detJ = np.linalg.det(J)
|
|
#print (f'detJ:{detJ}')
|
|
#print (J)
|
|
invJ = np.linalg.inv(J)
|
|
Jtrans=np.zeros((9+1,9+1))
|
|
Jtrans[1:4, 1:4]=invJ
|
|
Jtrans[4: 7, 4: 7]=invJ
|
|
Jtrans[7: 10, 7: 10]=invJ
|
|
Dif =np.zeros((9+1, 24+1))
|
|
for j in range(1,9):
|
|
for jj in range(1,4):
|
|
Dif[3 * jj - 2: 3 * jj+1, j * 3 + jj - 3] = dN[:, j-1]
|
|
|
|
|
|
B9x24 = np.matmul(Jtrans ,Dif)
|
|
B = np.zeros((6+1, 24+1))
|
|
B[1,:] = B9x24[1,:]
|
|
B[2,:] = B9x24[5,:]
|
|
B[3,:] = B9x24[9,:]
|
|
B[4,:] = B9x24[2,:] + B9x24[4,:]
|
|
B[5,:] = B9x24[6,:] + B9x24[8,:]
|
|
B[6,:] = B9x24[3,:]+ B9x24[7,:]
|
|
B=B[1:,1:]
|
|
|
|
deltaK=np.matmul(B.transpose(),C)
|
|
deltaK=np.matmul(deltaK,B)
|
|
K = K + deltaK*detJ
|
|
|
|
weight_x = (2 / (numgp_x * 2)) ** 3;
|
|
K = K*weight_x;
|
|
return K.reshape(24,24)
|
|
|
|
def test():
|
|
default = 0.123
|
|
x=[default]*10
|
|
|
|
DH=getDH(x)
|
|
print (DH)
|
|
print ('----------------')
|
|
X=np.asarray([0,0,0,0,1,1,1,1])
|
|
Y=np.asarray([0,0,1,1,0,0,1,1])
|
|
Z=np.asarray([0,1,0,1,0,1,0,1])
|
|
#X=np.random.rand(8)
|
|
#Y=np.random.rand(8)
|
|
#Z=np.random.rand(8)
|
|
KE=LocalSMRT(X,Y,Z,DH)
|
|
print (KE-KE.T)
|
|
|
|
|
|
|