pset 6

problem 1 part c

In [1]:
from sympy import *
In [2]:
h,x = symbols('h x')
M = Matrix([[1,0,0,0],[0,1,0,0],[1,h,h**2,h**3],[0,1,2*h,3*h**2]])
Minv = M.inv()
In [3]:
M
Out[3]:
$\displaystyle \left[\begin{matrix}1 & 0 & 0 & 0\\0 & 1 & 0 & 0\\1 & h & h^{2} & h^{3}\\0 & 1 & 2 h & 3 h^{2}\end{matrix}\right]$
In [4]:
Minv
Out[4]:
$\displaystyle \left[\begin{matrix}1 & 0 & 0 & 0\\0 & 1 & 0 & 0\\- \frac{3}{h^{2}} & - \frac{2}{h} & \frac{3}{h^{2}} & - \frac{1}{h}\\\frac{2}{h^{3}} & \frac{1}{h^{2}} & - \frac{2}{h^{3}} & \frac{1}{h^{2}}\end{matrix}\right]$
In [5]:
C = Matrix([1,x,x**2,x**3]) #entries for order three poly

u = Minv.T * C
du = diff(Minv.T * C,x)
In [6]:
u
Out[6]:
$\displaystyle \left[\begin{matrix}1 - \frac{3 x^{2}}{h^{2}} + \frac{2 x^{3}}{h^{3}}\\x - \frac{2 x^{2}}{h} + \frac{x^{3}}{h^{2}}\\\frac{3 x^{2}}{h^{2}} - \frac{2 x^{3}}{h^{3}}\\- \frac{x^{2}}{h} + \frac{x^{3}}{h^{2}}\end{matrix}\right]$
In [7]:
du
Out[7]:
$\displaystyle \left[\begin{matrix}- \frac{6 x}{h^{2}} + \frac{6 x^{2}}{h^{3}}\\1 - \frac{4 x}{h} + \frac{3 x^{2}}{h^{2}}\\\frac{6 x}{h^{2}} - \frac{6 x^{2}}{h^{3}}\\- \frac{2 x}{h} + \frac{3 x^{2}}{h^{2}}\end{matrix}\right]$

u basis

In [8]:
Aii = integrate(u[2]*u[2],(x,0,h)) + integrate(u[0]*u[0],(x,0,h)) #diagonal
Aii
Out[8]:
$\displaystyle \frac{26 h}{35}$
In [9]:
Aii_1  = integrate(u[0]*u[2],(x,0,h)) #u_i-1, u_i and the corresponding
Aii_1
Out[9]:
$\displaystyle \frac{9 h}{70}$
In [10]:
Bii = integrate(du[2]*du[2],(x,0,h)) + integrate(du[0]*du[0],(x,0,h)) #diagonal
Bii
Out[10]:
$\displaystyle \frac{12}{5 h}$
In [11]:
Bii_1 = integrate(du[0]*du[2],(x,0,h)) #du_i-1, du_i and the corresponding
Bii_1
Out[11]:
$\displaystyle - \frac{6}{5 h}$

du basis

In [12]:
Aii = integrate(u[3]*u[3],(x,0,h)) + integrate(u[1]*u[1],(x,0,h)) #diagonal
Aii
Out[12]:
$\displaystyle \frac{2 h^{3}}{105}$
In [13]:
Aii_1  = integrate(u[1]*u[3],(x,0,h)) #u_i-1, u_i and the corresponding
Aii_1
Out[13]:
$\displaystyle - \frac{h^{3}}{140}$
In [14]:
Bii = integrate(du[3]*du[3],(x,0,h)) + integrate(du[1]*du[1],(x,0,h)) #diagonal
Bii
Out[14]:
$\displaystyle \frac{4 h}{15}$
In [15]:
Bii_1 = integrate(du[1]*du[3],(x,0,h)) #du_i-1, du_i and the corresponding
Bii_1
Out[15]:
$\displaystyle - \frac{h}{30}$

u du

In [16]:
Aii = integrate(u[2]*u[3],(x,0,h)) + integrate(u[0]*u[1],(x,0,h)) #diagonal
Aii
Out[16]:
$\displaystyle 0$
In [17]:
Aii_1  = integrate(u[2]*u[1],(x,0,h)) #off diagonal
Aii_1
Out[17]:
$\displaystyle \frac{13 h^{2}}{420}$
In [18]:
Bii = integrate(du[2]*du[3],(x,0,h)) + integrate(du[0]*du[1],(x,0,h)) #diagonal
Bii
Out[18]:
$\displaystyle 0$
In [19]:
Bii_1 = integrate(du[1]*du[2],(x,0,h)) #du_i-1, du_i and the corresponding
Bii_1
Out[19]:
$\displaystyle - \frac{1}{10}$

problem 2

beam bending

In [2]:
import numpy as np
import matplotlib.pyplot as plt
In [3]:
#we reuse M, Minv, u, du from the the prior problem: replicated below:

h,x = symbols('h x')

M = Matrix([[1,0,0,0],[0,1,0,0],[1,h,h**2,h**3],[0,1,2*h,3*h**2]])
Minv = M.inv()

C = Matrix([1,x,x**2,x**3]) #entries for order three poly

u = Minv.T * C
du = diff(u, x)
ddu = diff(du, x)

#element stiffness matrix
k = zeros(4,4)

for i in range(k.shape[0]):
    for j in range(k.shape[1]):
        k[i,j] = integrate(ddu[i]*ddu[j],(x,0,h)).subs(h,1)
In [4]:
# node amount
N = 20
In [5]:
# below logic updated from Neil's example @ https://gitlab.cba.mit.edu/classes/864.23/site/-/blob/main/code/fem/beam.py
# bc i couldn't get it to work alas!
K = zeros(N,N)
K[0,0] = k[0,0]
K[0,1] = k[0,1]
K[0,2] = k[0,2]
K[0,3] = k[0,3]
K[1,0] = k[1,0]
K[1,1] = k[1,1]
K[1,2] = k[1,2]
K[1,3] = k[1,3]
for i in range(2,N-3,2):
    K[i,i-2] = k[2,0]
    K[i,i-1] = k[2,1]
    K[i,i] = k[2,2]+k[0,0]
    K[i,i+1] = k[2,3]+k[0,1]
    K[i,i+2] = k[0,2]
    K[i,i+3] = k[0,3]
    K[i+1,i-2] = k[3,0]
    K[i+1,i-1] = k[3,1]
    K[i+1,i] = k[3,2]+k[1,0]
    K[i+1,i+1] = k[3,3]+k[1,1]
    K[i+1,i+2] = k[1,2]
    K[i+1,i+3] = k[1,3]
K[N-2,N-4] = k[2,0]
K[N-2,N-3] = k[2,1]
K[N-2,N-2] = k[2,2]
K[N-2,N-1] = k[2,3]
K[N-1,N-4] = k[3,0]
K[N-1,N-3] = k[3,1]
K[N-1,N-2] = k[3,2]
K[N-1,N-1] = k[3,3]

K = np.array(K).astype(np.float64)

# boundary conditions
K[:,0] = 0
K[0,:] = 0
K[:,1] = 0
K[1,:] = 0
K[0,0] = 1
K[1,1] = 1

Kinv = np.linalg.inv(K)
In [6]:
from matplotlib.widgets import Slider
In [7]:
%matplotlib notebook
f_init = -1 # initial value
F = np.zeros(N)
F[N-1] = f_init

def f(force):
    F = np.zeros(N)
    F[N-1] = force
    return np.dot(Kinv,F)[::2]

fig, ax = plt.subplots()
line, = ax.plot(f(f_init))

# adjust the main plot to make room for the sliders
fig.subplots_adjust(left=0.25, bottom=0.25)

axForce = fig.add_axes([0.1, 0.25, 0.0225, 0.63])
force_Slider = Slider(
    ax=axForce,
    label="applied force",
    valmin=-10,
    valmax=0,
    valinit=f_init,
    orientation="vertical"
)

# The function to be called anytime a slider's value changes
def update(val):
    line.set_ydata(f(force_Slider.val))
    fig.canvas.draw_idle()

force_Slider.on_changed(update)

plt.show()
In [26]:
In [ ]:
# import taichi as ti
# ti.init(arch=ti.cpu)

# #construct global stiffness matrix

# K = ti.field(float, shape=(N,N))
# K.fill(0)

# @ti.kernel
# def fillK():
#     for i,j in K:
#         if i <=1 and j <= 3:
#             K[i,j] = k[i,j]
#         elif i >= (N-2) and j <=3:
#             K[i,(N-4+j)] = k[-(i-N),j]
#         elif j == i:
#             K[i,j] = k[0,0] + k[2,2]
#         else:
#             #ahugh
# fillK()
In [ ]:
import taichi as ti
ti.init(arch=ti.cpu)

# material properties
E = 2e9 # Pa
I = 100 # m^4

nodes = 10 # node amount
l = 0.1 # element length

p = ti.Vector.field(2, dtype=float, shape=(nodes,)) # x,y position of each node
theta = ti.Vector.field(2, dtype=float, shape=(nodes,)) # slope at each node

Fs = ti.Vector.field(1, dtype=float, shape=(nodes,)) # reaction shear at each node
M = ti.Vector.field(1, dtype=float, shape=(nodes,)) # moment at each node

u = ti.Vector.field(1, dtype=float, shape=(nodes,)) # applied force

# element stiffness matrix
kn = (E*I/l**2)*np.array([[12, 6*l, -12, 6*l],
        [6*l, 4*l**2, -6*l, 2*l**2],
        [-12, -6*l, 12, -6*l],
        [6*l, 2*l**2, -6*l, 4*l**4]])
k = ti.field(float, shape=(4,4))
k.from_numpy(kn)

# global stiffness matrix
K = ti.field(float, shape=(4*nodes,4*nodes))

@ti.kernel
def initBeam():
    for i in p:
        p[i] = [i*l, 0]
        theta = [i*l, 0]

@ti.kernel
def initU():
    for i in u:
        u[i] = [0]

@ti.kernel
def getGlobalStiffness():
    for i,j in K:
        K[i,j] #will finish later. 


#print(k)
initBeam()
initU()