deform.py

deform.py

import numpy as np
np.set_printoptions(precision=3)

# inputs
N = 2
E = 100
S = 0.5
nu = 0.2
A = np.array([[1,2,4],[4,3,1]])
xy = np.array([[0,0],[1,0],[0,1],[1,1]]) # coordinates [xi,yi]
Fe = [np.zeros([6]), np.array([0,-0.5,0,-0.5,0,0])]

nLocal = 3
nGlobal = 4

# D
D = np.array([[1-nu,nu,0],[nu,1-nu,0],[0,0,(1-2*nu)/2]])
print("#"*80)
print("Outputs:")
print("#"*80)
print("D:")
print(D)
print("-"*80)

# construct Pe for all elements
A -= 1 # python index from 0
Pe = []
for n in range(N):
    Pe_n = np.zeros([nLocal*2,nGlobal*2])
    for l,g in enumerate(A[n,:]):
        print("element {}: local {} <-> global {}".format(n+1,l+1,g+1))
        Pe_n[l*2:(l+1)*2, g*2:(g+1)*2] = np.eye(2)
    Pe.append(Pe_n)
    print("Pe_{}:".format(n+1))
    print(Pe[n].astype(np.int32))
print("-"*80)

# construct B_n from partial derivatives of N_n(x,y)
B = []
K = []
for n in range(N):
    xy_n = xy[A[n,:],:]
    xi = xy_n[0,0]
    xj = xy_n[1,0]
    xm = xy_n[2,0]
    yi = xy_n[0,1]
    yj = xy_n[1,1]
    ym = xy_n[2,1]
    delta = (xj-xi)*(ym-yi) - (yj-yi)*(xm-xi)
    dxN1 = -1/delta * (ym-yj)
    dxN2 =  1/delta * (ym-yi)
    dxN3 = -1/delta * (yj-yi)
    dyN1 = -1/delta * (xj-xm)
    dyN2 = -1/delta * (xm-xi)
    dyN3 =  1/delta * (xj-xi)
    B_n = np.array([[dxN1,0,dxN2,0,dxN3,0], [0,dyN1,0,dyN2,0,dyN3], [dyN1,dxN1,dyN2,dxN2,dyN3,dxN3]])
    DB = np.dot(D,B_n)
    Re = np.dot(B_n.transpose(), DB)
    K_n = np.dot(Pe[n].transpose(), np.dot(Re, Pe[n])) * E*S / ((1+nu)*(1-2*nu))
    B.append(B_n)
    K.append(K_n)
    print("B_{}:".format(n+1))
    print(B[n])
    print("Re_{}:".format(n+1))
    print(Re)
    print("K_{}:".format(n+1))
    print(K[n])
print("-"*80)

# RG
RG = np.zeros_like(K[0])
for n in range(N):
    RG += K[n]
print("RG:")
print(RG)
print("-"*80)

# F
F = np.zeros([Pe[0].shape[1]])
for n in range(N):
    F += np.dot(Pe[n].transpose(), Fe[n])
print("F:")
print(F)
print("-"*80)

# u
non_zeros = F != 0
rg = RG[:,non_zeros][non_zeros,:]
fnz = F[non_zeros]
u_nz = np.linalg.solve(rg, fnz)
u = np.zeros_like(F)
u[non_zeros] = u_nz
print("u:")
print(u)
print("-"*80)

# displacement
Eps = []
for n in range(N):
    eps_n = np.dot(B[n], np.dot(Pe[n], u))
    Eps.append(eps_n)
    print("Eps_{}:".format(n+1))
    print(Eps[n])
print("-"*80)

# deformation
G = []
for n in range(N):
    G_n = E / ((1+nu)*(1-2*nu)) * np.dot(D, Eps[n])
    G.append(G_n)
    print("G_{}:".format(n+1))
    print(G[n])
print("#"*80)