Implementation of dynamic system functions and remaining MATLAB examples

examples/exd_beam2_b.py

@@ -0,0 +1,80 @@
+# example exd_beam2_b
+# ----------------------------------------------------------------
+# PURPOSE
+#    Structural Dynamics, time integration, time dependent
+#    boundary conditions.
+#
+#    Note: file exd_beam2_m.py must be in the same directory
+# ----------------------------------------------------------------
+from exd_beam2_m import *
+
+# ----- Impact, center point, vertical beam ----------------------
+dt = 0.002
+T = 1
+
+# ----- Boundary condition, initial condition --------------------
+G = np.array([
+    [0, 0],
+    [0.1, 0.02],
+    [0.2, -0.01],
+    [0.3, 0],
+    [T, 0]
+])
+
+t,g = cfc.gfunc(G,dt)
+
+bc = np.zeros((4, 1+len(g)))
+bc[0,:] = np.hstack((1,g))
+bc[1,0] = 2
+bc[2,0] = 3
+bc[3,0] = 14
+
+a0 = np.zeros((15,1))
+da0 = np.zeros((15,1))
+
+# ----- Output parameters ----------------------------------------
+times = np.arange(0.1,1.1,0.1)
+dofs = np.array([1,4,11])
+
+# ----- Time integration parameters ------------------------------
+ip = np.array([dt, T, 0.25, 0.5])
+
+# ----- Time integration -----------------------------------------
+sol, dofhist = cfc.step2(K,[],M,[],a0,da0,bc,ip,times,dofs)
+
+# ----- Plot time history for two DOFs ---------------------------
+cfv.figure(1,fig_size=(7,4))
+cfv.plt.plot(t,dofhist['a'][0,:], '-')
+cfv.plt.plot(t,dofhist['a'][1,:],'--')
+cfv.plt.plot(t,dofhist['a'][2,:],'-.')
+cfv.plt.xlim([0,1])
+cfv.plt.ylim([-0.02,0.03])
+cfv.plt.xlabel("time (sec)")
+cfv.plt.ylabel("displacement (m)")
+cfv.plt.title("Displacement(time) at the 1st, 4th and 11th degree of freedom")
+cfv.text("solid line = bottom, vertical beam, x-direction", [0.2,0.022])
+cfv.text("dashed line = center, vertical beam, x-direction", [0.2,0.017])
+cfv.text("dashed-dotted line = center, horizontal beam, y-direction", [0.2,0.012])
+cfv.plt.grid()
+
+# ----- Plot displacement for some time increments ----------------
+cfv.figure(2,fig_size=(7,5))
+for i in range(5):
+    Edb = cfc.extract_ed(edof,sol['a'][:,i])
+    ext = ex+i*3.5
+    cfv.eldraw2(ext,ey,[2,3,1])
+    cfv.eldisp2(ext,ey,Edb,[1,2,2],sfac=20)
+    cfv.text(f"{times[i]:.1f}", [3.5*i+0.5,1.5])
+eyt = ey-4
+for i in range(5,10):
+    Edb = cfc.extract_ed(edof,sol['a'][:,i])
+    ext = ex+(i-5)*3.5
+    cfv.eldraw2(ext,eyt,[2,3,1])
+    cfv.eldisp2(ext,eyt,Edb,[1,2,2],sfac=20)
+    cfv.text(f"{times[i]:.1f}", [3.5*(i-5)+0.5,-2.5])
+cfv.title("Snapshots (sec), magnification = 20")
+ax = cfv.gca()
+ax.set_axis_off()
+cfv.showAndWait()
+
+# ----- End -------------------------------------------------------

examples/exd_beam2_m.py

@@ -0,0 +1,104 @@
+# example exd_beam2_m
+# ----------------------------------------------------------------
+# PURPOSE
+#    Set up the fe-model and perform eigenvalue analysis
+#    for a simple frame structure.
+# ----------------------------------------------------------------
+
+
+import numpy as np
+import calfem.core as cfc
+import calfem.vis_mpl as cfv
+
+# ----- Generate the model ---------------------------------------
+# ----- Material data --------------------------------------------
+E = 3e10
+Av=0.1030e-2
+Ah=0.0764e-2
+rho=2500
+Iv = 0.0171e-4
+Ih=0.00801e-4
+ep1 = [E, Av, Iv, rho*Av]   #IPE100
+ep2 = [E, Ah, Ih, rho*Ah]   #IPE80
+
+# ----- Topology -------------------------------------------------
+edof = np.array([
+    [1,2,3,4,5,6],
+    [4,5,6,7,8,9],
+    [7,8,9,10,11,12],
+    [10,11,12,13,14,15]
+])
+
+# ----- List of coordinates --------------------------------------
+coord = np.array([
+    [0.0, 0.0],
+    [0.0, 1.5],
+    [0.0, 3.0],
+    [1.0, 3.0],
+    [2.0, 3.0],
+])
+
+# ----- List of degrees of freedom -------------------------------
+dof = np.array([
+    [1, 2, 3],
+    [4, 5, 6],
+    [7, 8, 9],
+    [10, 11, 12],
+    [13, 14, 15],
+])
+
+# ----- Generate element matrices, assemble in global matrices ---
+K = np.zeros([15,15])
+M = np.zeros([15,15])
+
+ex, ey = cfc.coordxtr(edof, coord, dof);
+ep = np.array([ep1, ep1, ep2, ep2])
+
+for elx, ely, eltopo, elprop in zip(ex, ey, edof, ep):
+    Ke, Me = cfc.beam2d(elx, ely, elprop)
+    cfc.assem(eltopo, K, Ke)
+    cfc.assem(eltopo, M, Me)
+
+# ----- Eigenvalue analysis --------------------------------------
+b = np.array([1, 2, 3, 14])
+La, Egv = cfc.eigen(K,M,b)
+freq = np.sqrt(La)/(2*np.pi)
+
+if __name__=='__main__':
+    # ----- Draw a plot of the element mesh --------------------------
+    cfv.figure(1,fig_size=(5.5, 4.5))
+    cfv.eldraw2(ex,ey,[1, 2, 1])
+    cfv.title('2-D Frame Structure')
+
+    # ----- Plot one eigenmode ---------------------------------------
+    cfv.figure(2,fig_size=(5.5, 4.5))
+    cfv.eldraw2(ex,ey,[2, 3, 1])
+    Edb = cfc.extract_ed(edof, Egv[:,0])
+    cfv.eldisp2(ex,ey,Edb,[1, 2, 2])
+    cfv.title('The first eigenmode')
+    cfv.text(f"{freq[0]:.2f}", [0.5,1.75])
+    ax = cfv.gca()
+    ax.grid()
+
+    # ----- Plot eight eigenmodes ------------------------------------
+    cfv.figure(3,fig_size=(7,5))
+    for i in range(4):
+        Edb = cfc.extract_ed(edof,Egv[:,i])
+        ext = ex+i*3
+        cfv.eldraw2(ext,ey,[2,3,1])
+        cfv.eldisp2(ext,ey,Edb,[1,2,2],sfac=0.5)
+        cfv.text(f"{freq[i]:.2f}", [3*i+0.5,1.5])
+    eyt = ey-4
+    for i in range(4,8):
+        Edb = cfc.extract_ed(edof,Egv[:,i])
+        ext = ex+(i-4)*3
+        cfv.eldraw2(ext,eyt,[2,3,1])
+        cfv.eldisp2(ext,eyt,Edb,[1,2,2],sfac=0.5)
+        cfv.text(f"{freq[i]:.2f}", [3*(i-4)+0.5,-2.5])
+    cfv.title("The first eight eigenmodes [Hz]")
+    ax = cfv.gca()
+    ax.set_axis_off()
+    cfv.showAndWait()
+
+# ----- End -------------------------------------------------------
+

examples/exd_beam2_t.py

@@ -0,0 +1,80 @@
+# example exd_beam2_t
+# ----------------------------------------------------------------
+# PURPOSE
+#    Structural Dynamics, time integration, full system.
+#
+#    Note: file exd_beam2_m.py must be in the same directory
+# ----------------------------------------------------------------
+from exd_beam2_m import *
+
+# ----- Impact, center point, vertical beam ----------------------
+dt = 0.002
+T = 1
+
+# ----- The load -------------------------------------------------
+G = np.array([
+    [0, 0],
+    [0.15, 1],
+    [0.25, 0],
+    [T, 0]
+])
+
+t,g = cfc.gfunc(G,dt)
+f = np.zeros((15, len(g)))
+f[3,:] = 1000*g
+
+# ----- Boundary condition, initial condition --------------------
+bc = np.array([
+    [1, 0],
+    [2, 0],
+    [3, 0],
+    [14, 0]
+])
+
+a0 = np.zeros((15,1))
+da0 = np.zeros((15,1))
+
+# ----- Output parameters ----------------------------------------
+times = np.arange(0.1,1.1,0.1)
+dofs = np.array([4,11])
+
+# ----- Time integration parameters ------------------------------
+ip = np.array([dt, T, 0.25, 0.5])
+
+# ----- Time integration -----------------------------------------
+sol, dofhist = cfc.step2(K,[],M,f,a0,da0,bc,ip,times,dofs)
+
+# ----- Plot time history for two DOFs ---------------------------
+cfv.figure(1,fig_size=(7,4))
+cfv.plt.plot(t,dofhist['a'][0,:], '-')
+cfv.plt.plot(t,dofhist['a'][1,:],'--')
+cfv.plt.xlim([0,1])
+cfv.plt.ylim([-0.01,0.02])
+cfv.plt.xlabel("time (sec)")
+cfv.plt.ylabel("displacement (m)")
+cfv.plt.title("Displacement(time) at the 4th and 11th degree of freedom")
+cfv.text("solid line = impact point, x-direction", [0.3,0.017])
+cfv.text("dashed line = center, horizontal beam, y-direction", [0.3,0.012])
+cfv.plt.grid()
+
+# ----- Plot displacement for some time increments ----------------
+cfv.figure(2,fig_size=(7,5))
+for i in range(5):
+    Edb = cfc.extract_ed(edof,sol['a'][:,i])
+    ext = ex+i*3.5
+    cfv.eldraw2(ext,ey,[2,3,1])
+    cfv.eldisp2(ext,ey,Edb,[1,2,2],sfac=25)
+    cfv.text(f"{times[i]:.1f}", [3.5*i+0.5,1.5])
+eyt = ey-4
+for i in range(5,10):
+    Edb = cfc.extract_ed(edof,sol['a'][:,i])
+    ext = ex+(i-5)*3.5
+    cfv.eldraw2(ext,eyt,[2,3,1])
+    cfv.eldisp2(ext,eyt,Edb,[1,2,2],sfac=25)
+    cfv.text(f"{times[i]:.1f}", [3.5*(i-5)+0.5,-2.5])
+cfv.title("Snapshots (sec), magnification = 25")
+ax = cfv.gca()
+ax.set_axis_off()
+cfv.showAndWait()
+
+# ----- End -------------------------------------------------------

examples/exd_beam2_tr.py

@@ -0,0 +1,85 @@
+# example exd_beam2_tr
+# ----------------------------------------------------------------
+# PURPOSE
+#    Structural Dynamics, time integration, reduced system.
+#
+#    Note: file exd_beam2_m.py must be in the same directory
+# ----------------------------------------------------------------
+from exd_beam2_m import *
+
+# ----- Impact, center point, vertical beam ----------------------
+dt = 0.002
+T = 1
+nev = 2
+
+# ----- The load -------------------------------------------------
+G = np.array([
+    [0, 0],
+    [0.15, 1],
+    [0.25, 0],
+    [T, 0]
+])
+
+t,g = cfc.gfunc(G,dt)
+f = np.zeros((15, len(g)))
+f[3,:] = 1000*g
+fr = np.hstack((np.arange(1,nev+1).reshape(-1,1),Egv[:,:nev].T@f))
+
+# ----- Reduced system matrices ----------------------------------
+kr = np.diag(np.diag(Egv[:,:nev].T@K@Egv[:,:nev]))
+mr = np.diag(np.diag(Egv[:,:nev].T@M@Egv[:,:nev]))
+
+# ----- Initial condition ----------------------------------------
+ar0 = np.zeros((nev,1))
+dar0 = np.zeros((nev,1))
+
+# ----- Output parameters ----------------------------------------
+times = np.arange(0.1,1.1,0.1)
+dofsr = np.arange(1,nev+1)
+dofs = np.array([4,11])
+
+# ----- Time integration parameters ------------------------------
+ip = np.array([dt, T, 0.25, 0.5])
+
+# ----- Time integration -----------------------------------------
+sol, dofhist = cfc.step2(kr,[],mr,fr,ar0,dar0,[],ip,times,dofsr)
+
+# ----- Mapping back to original coordinate system ---------------
+aR = Egv[:,:nev]@sol['a']
+aRhist = Egv[dofs-1,:nev]@dofhist['a']
+
+# ----- Plot time history for two DOFs ---------------------------
+cfv.figure(1,fig_size=(7,4))
+cfv.plt.plot(t,aRhist[0,:], '-')
+cfv.plt.plot(t,aRhist[1,:],'--')
+cfv.plt.xlim([0,1])
+cfv.plt.ylim([-0.01,0.02])
+cfv.plt.xlabel("time (sec)")
+cfv.plt.ylabel("displacement (m)")
+cfv.plt.title("Displacement(time) at the 4th and 11th degree of freedom")
+cfv.text("solid line = impact point, x-direction", [0.3,0.017])
+cfv.text("dashed line = center, horizontal beam, y-direction", [0.3,0.012])
+cfv.text("TWO EIGENVECTORS ARE USED", [0.3,-0.007])
+cfv.plt.grid()
+
+# ----- Plot displacement for some time increments ----------------
+cfv.figure(2,fig_size=(7,5))
+for i in range(5):
+    Edb = cfc.extract_ed(edof,aR[:,i])
+    ext = ex+i*3.5
+    cfv.eldraw2(ext,ey,[2,3,1])
+    cfv.eldisp2(ext,ey,Edb,[1,2,2],sfac=25)
+    cfv.text(f"{times[i]:.1f}", [3.5*i+0.5,1.5])
+eyt = ey-4
+for i in range(5,10):
+    Edb = cfc.extract_ed(edof,aR[:,i])
+    ext = ex+(i-5)*3.5
+    cfv.eldraw2(ext,eyt,[2,3,1])
+    cfv.eldisp2(ext,eyt,Edb,[1,2,2],sfac=25)
+    cfv.text(f"{times[i]:.1f}", [3.5*(i-5)+0.5,-2.5])
+cfv.title("Snapshots (sec), magnification = 25")
+ax = cfv.gca()
+ax.set_axis_off()
+cfv.showAndWait()
+
+# ----- End -------------------------------------------------------