From 796587f85d1aa0f5962039da06171b8e7f82da5c Mon Sep 17 00:00:00 2001 From: Adam Sciegaj Date: Thu, 15 Dec 2022 12:49:07 +0100 Subject: [PATCH 1/7] Recreated example exd_beam2_m --- examples/exd_beam2_m.py | 104 ++++++++++++++++++++++++++++++++++++++++ 1 file changed, 104 insertions(+) create mode 100644 examples/exd_beam2_m.py diff --git a/examples/exd_beam2_m.py b/examples/exd_beam2_m.py new file mode 100644 index 0000000..266fccd --- /dev/null +++ b/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) + +# ----- 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') + +# ----- Eigenvalue analysis -------------------------------------- +b = np.array([1, 2, 3, 14]) +La, Egv = cfc.eigen(K,M,b) +freq = np.sqrt(La)/(2*np.pi) + +# ----- 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 ------------------------------------------------------- + From 16c72f92e174d2ab2c24343edc7ae7ffb4c54243 Mon Sep 17 00:00:00 2001 From: Adam Sciegaj Date: Fri, 16 Dec 2022 15:53:51 +0100 Subject: [PATCH 2/7] step2 and gfunc functions --- calfem/core.py | 228 ++++++++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 227 insertions(+), 1 deletion(-) diff --git a/calfem/core.py b/calfem/core.py index b4e58c0..ff70733 100644 --- a/calfem/core.py +++ b/calfem/core.py @@ -6,7 +6,7 @@ """ from scipy.sparse.linalg import dsolve -from scipy.linalg import eig +from scipy.linalg import eig, lu import numpy as np import logging as cflog @@ -3936,6 +3936,232 @@ def eigen(K,M,b=None): return L, X +def gfunc(G,dt): + """ + Form vector with function values at equally spaced + points by linear interpolation + + Parameters: + + G = [t_i, g_i] t_i: time i, g_i: g(t_i) + dim(G) = np x 2, np = number of points + dt time step + + Returns: + + t 1-D vector with equally spaced time points + g 1-D vector with corresponding function values + """ + ti = np.arange(G[0,0],G[-1,0]+dt,dt) + g1 = np.interp(ti,G[:,0],G[:,1]) + return ti, g1 + + +def step2(K,C,M,f,a0,da0,bc,ip,times,dofs): + """ + Algorithm for dynamic solution of second-order + FE equations considering boundary conditions. + + Parameters: + + K global stiffness matrix, dim(K) = ndof x ndof + C global damping matrix, dim(C) = ndof x ndof + If there is no damping in the system, simply set C=[] + M global mass matrix, dim(M) = ndof x ndof + f global load vector, dim(f) = ndof x (nstep + 1), + If dim(f) = ndof x 1, the values are kept constant + during time integration + a0 initial displacement vector a(0), dim(a0) = ndof x 1 + da0 initial velocity vector v(0), dim(da0) = ndof x 1 + bc boundary condition matrix, dim(bc) = nbc x (nstep + 2) + where nbc = number of prescribed degrees of freedom (either constant or time-dependent) + The first column contains the numbers of the prescribed degrees of freedom + and the subsequent columns contain the time history. + If dim(bc) = nbc x 2, the values from the second column are kept constant + during time integration + ip list [dt, tottime, alpha, delta], where + dt is the size of the time increment, + tottime is the total time, + alpha and delta are time integration constants for the Newmark family of methods. + Frequently used values of alpha and delta: + alpha=1/4, delta=1/2: average acceleration (trapezoidal) rule, + alpha=1/6, delta=1/2: linear acceleration + alpha=0, delta=1/2: central difference + times array [t(i) ...] of times at which output should be written to a, da and d2a + dofs array [dof(i) ...] of degree of freedom numbers for which history output + should be written to ahist, dahist and d2ahist + + Returns: + + modelhist dictionary containing solution history for the whole model at following keys: + modelhist['a'] constains displacement values at all timesteps, + alternatively at times specified in 'times' + dim(modelhist['a']) = ndof x (nstep + 1) or ndof x ntimes + modelhist['da'] constains velocity values at all timesteps, + alternatively at times specified in 'times' + dim(modelhist['da']) = ndof x (nstep + 1) or ndof x ntimes + modelhist['d2a'] constains acceleration values at all timesteps, + alternatively at times specified in 'times' + dim(modelhist['d2a']) = ndof x (nstep + 1) or ndof x ntimes + dofhist dictionary containing solution history for the selected 'dofs': + dofhist['a'] constains displacement time history at the dofs specified in 'dofs' + dim(dofhist['ahist']) = ndof x (nstep + 1) + dofhist['da'] constains velocity time history at the dofs specified in 'dofs' + dim(dofhist['dahist']) = ndof x (nstep + 1) + dofhist['d2a'] constains acceleration time history at the dofs specified in 'dofs' + dim(dofhist['d2ahist']) = ndof x (nstep + 1) + """ + ndof, _ = K.shape + if not C: + C = np.zeros((ndof,ndof)) + dt, tottime, alpha, delta = ip + b1 = dt*dt*0.5*(1-2*alpha) + b2 = (1-delta)*dt + b3 = delta*dt + b4 = alpha*dt*dt + + nstep = 1 + if np.array(f).any(): + _, ncf = f.shape + if ncf>1: + nstep = ncf-1 + + if np.array(bc).any(): + _, ncb = bc.shape + if ncb>2: + nstep = ncb-2 + bound = 1 + if not np.array(bc).any(): + bound = 0 + + ns = int(tottime/dt) + if (ns < nstep or nstep==1): + nstep=ns + + tf = np.zeros((ndof,nstep+1)) + if np.array(f).any(): + if ncf==1: + tf = f[:,0].reshape(-1,1)@np.ones((1,nstep+1)) + if ncf>1: + tf = np.copy(f) + + modelhist = {} + sa=0 + if not times.any(): + ntimes=0 + sa=1 + modelhist['a'] = np.zeros((ndof,nstep+1)) + modelhist['da'] = np.zeros((ndof,nstep+1)) + modelhist['d2a'] = np.zeros((ndof,nstep+1)) + else: + ntimes = len(times) + if ntimes: + sa=2 + modelhist['a'] = np.zeros((ndof,ntimes)) + modelhist['da'] = np.zeros((ndof,ntimes)) + modelhist['d2a'] = np.zeros((ndof,ntimes)) + + dofhist = {} + if dofs.all(): + ndofs = len(dofs) + if ndofs: + dofhist['a'] = np.zeros((ndofs,nstep+1)) + dofhist['da'] = np.zeros((ndofs,nstep+1)) + dofhist['d2a'] = np.zeros((ndofs,nstep+1)) + else: + ndofs=0 + + itime = 0 + + # Calculate initial second time derivative d2a0 + d2a0 = np.linalg.solve(M,tf[:,0].reshape(-1,1) - C@da0 - K@a0) + # Save initial values + if sa==1: + modelhist['a'][:,0] = a0.ravel() + modelhist['da'][:,0] = da0.ravel() + modelhist['d2a'][:,0] = d2a0.ravel() + elif sa==2: + if times[itime]==0: + modelhist['a'][:,itime] = a0.ravel() + modelhist['da'][:,itime] = da0.ravel() + modelhist['d2a'][:,itime] = d2a0.ravel() + itime += 1 + + if ndofs: + dofhist['a'][:,0] = a0[np.ix_(dofs-1)].ravel() + dofhist['da'][:,0] = da0[np.ix_(dofs-1)].ravel() + dofhist['d2a'][:,0] = d2a0[np.ix_(dofs-1)].ravel() + + # Reduce matrices due to bcs + tempa = np.zeros((ndof,1)) + tempda = np.zeros((ndof,1)) + tempd2a = np.zeros((ndof,1)) + fdof=np.arange(1,ndof+1).astype(int) + if bound: + nrb, ncb = bc.shape + if ncb==2: + pa = bc[:,1].reshape(-1,1)@np.ones((1,nstep+1)) + pda = np.zeros((nrb,nstep+1)) + elif ncb>2: + pa = np.copy(bc[:,1:]) + pda1 = (pa[:,1]-pa[:,0])/dt + pdarest = (pa[:,1:] - pa[:,0:-1])/dt + pda = np.hstack((pda1.reshape(-1,1),pdarest)) + pdof = np.copy(bc[:,0]).astype(int) + fdof = np.setdiff1d(fdof,pdof).astype(int) - 1 + pdof -= 1 #adjusting for indexing starting from 0 + Keff = M[np.ix_(fdof,fdof)] + b3*C[np.ix_(fdof,fdof)] +b4*K[np.ix_(fdof,fdof)] + else: + fdof -= 1 #adjusting for indexing starting from 0 + Keff = M + b3*C + b4*K + + L, U = lu(Keff,permute_l=True) + anew = a0[np.ix_(fdof)] + danew = da0[np.ix_(fdof)] + d2anew = d2a0[np.ix_(fdof)] + + # Iterate over time steps + for j in range(1,nstep+1): + time = dt*j + aold = np.copy(anew) + daold = np.copy(danew) + d2aold = np.copy(d2anew) + apred = aold + dt*daold + b1*d2aold + dapred = daold + b2*d2aold + if not bound: + reff = tf[:,j].reshape(-1,1) - C@dapred - K@apred + else: + pdeff = C[np.ix_(fdof,pdof)]@pda[:,j].reshape(-1,1) + K[np.ix_(fdof,pdof)]@pa[:,j].reshape(-1,1) + reff = tf[np.ix_(fdof),j].reshape(-1,1) - C[np.ix_(fdof,fdof)]@dapred - K[np.ix_(fdof,fdof)]@apred - pdeff + y = np.linalg.solve(L,reff) + d2anew = np.linalg.solve(U,y) + anew = apred + b4*d2anew + danew = dapred + b3*d2anew + # Save to modelhist and dofhist + if bound: + tempa[np.ix_(pdof)] = pa[:,j].reshape(-1,1) + tempda[np.ix_(pdof)] = pda[:,j].reshape(-1,1) + tempa[np.ix_(fdof)] = anew + tempda[np.ix_(fdof)] = danew + tempd2a[np.ix_(fdof)] = d2anew + if sa==1: + modelhist['a'][:,j] = tempa.ravel() + modelhist['da'][:,j] = tempda.ravel() + modelhist['d2a'][:,j] = tempd2a.ravel() + elif sa==2: + if ntimes and itime < ntimes: + if time >= times[itime]: + modelhist['a'][:,itime] = tempa.ravel() + modelhist['da'][:,itime] = tempda.ravel() + modelhist['d2a'][:,itime] = tempd2a.ravel() + itime += 1 + if ndofs: + dofhist['a'][:,j] = tempa[np.ix_(dofs-1)].ravel() + dofhist['da'][:,j] = tempda[np.ix_(dofs-1)].ravel() + dofhist['d2a'][:,j] = tempd2a[np.ix_(dofs-1)].ravel() + + return modelhist, dofhist + def extract_eldisp(edof, a): """ Extract element displacements from the global displacement From b2761b7997bb68ed9bf93f60707a66b0f03567ae Mon Sep 17 00:00:00 2001 From: Adam Sciegaj Date: Fri, 16 Dec 2022 15:54:45 +0100 Subject: [PATCH 3/7] Added examples exd_beam2_t, exd_beam2_tr, and exd_beam2_b --- examples/exd_beam2_b.py | 80 +++++++++++++++++++++++++++++++++++++ examples/exd_beam2_m.py | 72 +++++++++++++++++----------------- examples/exd_beam2_t.py | 80 +++++++++++++++++++++++++++++++++++++ examples/exd_beam2_tr.py | 85 ++++++++++++++++++++++++++++++++++++++++ 4 files changed, 281 insertions(+), 36 deletions(-) create mode 100644 examples/exd_beam2_b.py create mode 100644 examples/exd_beam2_t.py create mode 100644 examples/exd_beam2_tr.py diff --git a/examples/exd_beam2_b.py b/examples/exd_beam2_b.py new file mode 100644 index 0000000..71dac27 --- /dev/null +++ b/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 ------------------------------------------------------- \ No newline at end of file diff --git a/examples/exd_beam2_m.py b/examples/exd_beam2_m.py index 266fccd..99be247 100644 --- a/examples/exd_beam2_m.py +++ b/examples/exd_beam2_m.py @@ -11,7 +11,7 @@ import calfem.vis_mpl as cfv # ----- Generate the model --------------------------------------- -# ----- Material data --------------------------------------------- +# ----- Material data -------------------------------------------- E = 3e10 Av=0.1030e-2 Ah=0.0764e-2 @@ -59,46 +59,46 @@ cfc.assem(eltopo, K, Ke) cfc.assem(eltopo, M, Me) -# ----- 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') - # ----- Eigenvalue analysis -------------------------------------- b = np.array([1, 2, 3, 14]) La, Egv = cfc.eigen(K,M,b) freq = np.sqrt(La)/(2*np.pi) -# ----- 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() +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 ------------------------------------------------------- diff --git a/examples/exd_beam2_t.py b/examples/exd_beam2_t.py new file mode 100644 index 0000000..6bac640 --- /dev/null +++ b/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 ------------------------------------------------------- \ No newline at end of file diff --git a/examples/exd_beam2_tr.py b/examples/exd_beam2_tr.py new file mode 100644 index 0000000..0dd929f --- /dev/null +++ b/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 ------------------------------------------------------- \ No newline at end of file From f9fc27f9ebf09bdd5bd6c6824f1a333f3c953767 Mon Sep 17 00:00:00 2001 From: Adam Sciegaj Date: Sat, 17 Dec 2022 18:25:02 +0100 Subject: [PATCH 4/7] Fixed a bug in beam1s. Added example exs_beam2 and plots in exs_bar2 and exs_beam1 --- calfem/core.py | 12 ++--- examples/exs_bar2.py | 26 ++++++++++ examples/exs_beam1.py | 85 ++++++++++++++++--------------- examples/exs_beam2.py | 114 ++++++++++++++++++++++++++++++++++++++++++ 4 files changed, 192 insertions(+), 45 deletions(-) create mode 100644 examples/exs_beam2.py diff --git a/calfem/core.py b/calfem/core.py index ff70733..ffa56c7 100644 --- a/calfem/core.py +++ b/calfem/core.py @@ -307,7 +307,7 @@ def beam1e(ex, ep, eq=None): :param list ep: element properties [E, I], E - Young's modulus, I - Moment of inertia :param float eq: distributed load [qy] :return mat Ke: element stiffness matrix [4 x 4] - :return mat fe: element stiffness matrix [4 x 1] (if eq!=None) + :return mat fe: element load vector [4 x 1] (if eq!=None) """ L = ex[1]-ex[0] @@ -315,7 +315,7 @@ def beam1e(ex, ep, eq=None): I = ep[1] qy = 0. - if not eq is None: + if eq: qy = eq Ke = E*I/(L**3) * np.mat([ @@ -370,7 +370,7 @@ def beam1s(ex, ep, ed, eq=None, nep=None): qy = 0. - if not eq is None: + if eq: qy = eq ne = 2 @@ -387,15 +387,15 @@ def beam1s(ex, ep, ed, eq=None, nep=None): Ca = (Cinv@ed).T - x = np.asmatrix(np.arange(0., L+L/(ne-1), L/(ne-1))).T + x = np.asmatrix(np.linspace(0., L, nep)).T zero = np.asmatrix(np.zeros([len(x)])).T one = np.asmatrix(np.ones([len(x)])).T v = np.concatenate((one, x, np.power(x, 2), np.power(x, 3)), 1)@Ca \ + qy/(24*EI)*(np.power(x,4) - 2*L*np.power(x,3) + (L**2)*np.power(x,2)) d2v = np.concatenate((zero, zero, 2*one, 6*x), 1)@Ca \ - + qy/(2*EI)*(np.power(x,2) - L*x + L**2/12) - d3v = np.concatenate((zero, zero, zero, 6*one), 1)@Ca - qy*(x - L/2) + + qy/(12*EI)*(6*np.power(x,2) - 6*L*x + L**2) + d3v = np.concatenate((zero, zero, zero, 6*one), 1)@Ca + qy/(2*EI)*(2*x - L) M = EI*d2v V = -EI*d3v diff --git a/examples/exs_bar2.py b/examples/exs_bar2.py index e332348..dbf9cbd 100644 --- a/examples/exs_bar2.py +++ b/examples/exs_bar2.py @@ -13,6 +13,7 @@ import numpy as np import calfem.core as cfc +import calfem.vis_mpl as cfv # ----- Topology matrix Edof ------------------------------------- @@ -86,3 +87,28 @@ print("N1 = "+str(N1)) print("N2 = "+str(N2)) print("N3 = "+str(N3)) + +# ----- Draw deformed truss -------------------------------------- +cfv.figure(1,fig_size=(6,4)) +plotpar = [2,1,0] +cfv.eldraw2(ex1,ey1,plotpar) +cfv.eldraw2(ex2,ey2,plotpar) +cfv.eldraw2(ex3,ey3,plotpar) +sfac = cfv.scalfact2(ex1,ey1,ed1,0.1) +plotpar = [1,2,1] +cfv.eldisp2(ex1,ey1,ed1,plotpar,sfac) +cfv.eldisp2(ex2,ey2,ed2,plotpar,sfac) +cfv.eldisp2(ex3,ey3,ed3,plotpar,sfac) +cfv.plt.axis([-0.4, 2.0, -0.4, 1.4]) +cfv.title("Displacements") + +# ----- Draw normal force diagram --------------------------------- +cfv.figure(2,fig_size=(6,4)) +plotpar = [2,1] +sfac = cfv.scalfact2(ex1,ey1,N2,0.1) +cfv.secforce2(ex1,ey1,N1,plotpar,sfac) +cfv.secforce2(ex2,ey2,N2,plotpar,sfac) +cfv.secforce2(ex3,ey3,N3,plotpar,sfac) +cfv.plt.axis([-0.4, 2.0, -0.4, 1.4]) +cfv.title("Normal force") +cfv.showAndWait() \ No newline at end of file diff --git a/examples/exs_beam1.py b/examples/exs_beam1.py index 2c557b1..fafe72d 100644 --- a/examples/exs_beam1.py +++ b/examples/exs_beam1.py @@ -14,67 +14,74 @@ import numpy as np import calfem.core as cfc +import calfem.vis_mpl as cfv # ----- Topology ------------------------------------------------- edof = np.array([ - [1, 2, 3, 4, 5, 6], - [4, 5, 6, 7, 8, 9], - [7, 8, 9, 10, 11, 12] + [1, 2, 3, 4], + [3, 4, 5, 6], ]) # ----- Stiffness matrix K and load vector f --------------------- -K = np.mat(np.zeros((12, 12))) -f = np.mat(np.zeros((12, 1))) -f[4] = -10000. +K = np.mat(np.zeros((6, 6))) +f = np.mat(np.zeros((6, 1))) +f[2] = -10000. # ----- Element stiffness matrices ------------------------------ E = 2.1e11 -A = 45.3e-4 I = 2510e-8 -ep = np.array([E, A, I]) -ex = np.array([0., 3.]) -ey = np.array([0., 0.]) - -Ke = cfc.beam2e(ex, ey, ep) - -print(Ke) +ep = np.array([E, I]) +ex1 = np.array([0., 3.]) +ex2 = np.array([3., 9.]) +Ke1 = cfc.beam1e(ex1, ep) +Ke2 = cfc.beam1e(ex2, ep) # ----- Assemble Ke into K --------------------------------------- -K = cfc.assem(edof, K, Ke) +K = cfc.assem(edof[0,:], K, Ke1) +K = cfc.assem(edof[1,:], K, Ke2) # ----- Solve the system of equations and compute support forces - -bc = np.array([1, 2, 11]) +bc = np.array([1, 5]) (a, r) = cfc.solveq(K, f, bc) # ----- Section forces ------------------------------------------- ed = cfc.extract_ed(edof, a) -es1, ed1, ec1 = cfc.beam2s(ex, ey, ep, ed[0, :], nep=10) -es2, ed2, ec2 = cfc.beam2s(ex, ey, ep, ed[1, :], nep=10) -es3, ed3, ec3 = cfc.beam2s(ex, ey, ep, ed[2, :], nep=10) - -# ----- Results -------------------------------------------------- - -print("a=") -print(a) -print("r=") -print(r) -print("es1=") -print(es1) -print("es2=") -print(es2) -print("es3=") -print(es3) - -print("ed1=") -print(ed1) -print("ed2=") -print(ed2) -print("ed3=") -print(ed3) +es1, ed1, ec1 = cfc.beam1s(ex1, ep, ed[0, :], nep=10) +es2, ed2, ec2 = cfc.beam1s(ex2, ep, ed[1, :], nep=10) + +es = np.vstack(([0,0],es1,es2,[0,0])) +eD = np.vstack((0,ed1,ed2,0)) +ec = np.vstack((ex1[0],ec1+ex1[0],ec2+ex2[0],ex2[1])) + +# ----- Draw deformed beam ---------------------------------------- +cfv.figure(1,fig_size=(6,2.5)) +cfv.plt.plot([0, 9],[0,0],'b',linewidth=0.5) +cfv.plt.plot(ec,eD,'b') +cfv.plt.axis([-1,10,-0.03, 0.01]) +cfv.title("Displacements") + +# ----- Draw shear force diagram ---------------------------------- +cfv.figure(2,fig_size=(6,2.5)) +cfv.plt.plot([0, 9],[0,0],'b',linewidth=0.5) +cfv.plt.plot(ec,es[:,0],'b') +cfv.plt.axis([-1,10,-8000, 5000]) +cfv.title("Shear force") +cfv.gca().invert_yaxis() + +# ----- Draw bending moment diagram ------------------------------- +cfv.figure(3,fig_size=(6,2.5)) +cfv.plt.plot([0, 9],[0,0],'b',linewidth=0.5) +cfv.plt.plot(ec,es[:,1],'b') +cfv.plt.axis([-1,10,-5000, 25000]) +cfv.title("Bending moment") +cfv.gca().invert_yaxis() +cfv.showAndWait() + +#------------------------ end ----------------------------------- \ No newline at end of file diff --git a/examples/exs_beam2.py b/examples/exs_beam2.py new file mode 100644 index 0000000..98869ac --- /dev/null +++ b/examples/exs_beam2.py @@ -0,0 +1,114 @@ +# example exs6 +# ---------------------------------------------------------------- +# PURPOSE +# Analysis of a plane frame. +# ---------------------------------------------------------------- + +import numpy as np +import calfem.core as cfc +import calfem.vis_mpl as cfv + + +#----- Topology ------------------------------------------------- + +edof = np.array([ + [4, 5, 6, 1, 2, 3], + [7, 8, 9, 10, 11, 12], + [4, 5, 6, 7, 8, 9] +]) + +#----- Stiffness matrix K and load vector f --------------------- + +K = np.zeros([12,12]) +f = np.zeros([12,1]) +f[3,0] = 2e3 + +#----- Element stiffness and element load matrices ------------- + +E=200e9 +A1=2e-3 +A2=6e-3 +I1=1.6e-5 +I2=5.4e-5 + +ep1 = np.array([E, A1, I1]) +ep3 = np.array([E, A2, I2]) + +ex1 = np.array([0.,0.]) +ey1 = np.array([4.,0.]) +ex2 = np.array([6.,6.]) +ey2 = np.array([4.,0.]) +ex3 = np.array([0.,6.]) +ey3 = np.array([4.,4.]) + +eq1 = np.array([0.,0.]) +eq2 = np.array([0.,0.]) +eq3 = np.array([0.,-10e3]) + +Ke1 = cfc.beam2e(ex1,ey1,ep1) +Ke2 = cfc.beam2e(ex2,ey2,ep1) +Ke3, fe3 = cfc.beam2e(ex3,ey3,ep3,eq3) + +#----- Assemble Ke into K --------------------------------------- + +K = cfc.assem(edof[0,:],K,Ke1) +K = cfc.assem(edof[1,:],K,Ke2) +K, f = cfc.assem(edof[2,:],K,Ke3,f,fe3) + +#----- Solve the system of equations and compute reactions ------ +bc = np.array([1, 2, 3, 10, 11]) +a, r = cfc.solveq(K, f, bc) + +#----- Section forces ------------------------------------------- +ed = cfc.extract_ed(edof, a) + +es1, edi1, eci1 = cfc.beam2s(ex1,ey1,ep1,ed[0,:],eq1,21) +es2, edi2, eci2 = cfc.beam2s(ex2,ey2,ep1,ed[1,:],eq2,21) +es3, edi3, eci3 = cfc.beam2s(ex3,ey3,ep3,ed[2,:],eq3,21) + +#----- Draw deformed frame --------------------------------------- +cfv.figure(1,fig_size=(6,4)) +plotpar=[2,1,0] +cfv.eldraw2(ex1,ey1,plotpar) +cfv.eldraw2(ex2,ey2,plotpar) +cfv.eldraw2(ex3,ey3,plotpar) +sfac=cfv.scalfact2(ex3,ey3,edi3,0.1) +plotpar=[1,2,1] +cfv.eldisp2(ex1,ey1,ed[0,:],plotpar,sfac) +cfv.eldisp2(ex2,ey2,ed[1,:],plotpar,sfac) +cfv.eldisp2(ex3,ey3,ed[2,:],plotpar,sfac) +cfv.plt.axis([-1.5, 7.5, -0.5, 5.5]) +cfv.title("Displacements") + + +#----- Draw normal force diagram -------------------------------- +cfv.figure(2,fig_size=(6,4)) +plotpar=[2,1] +sfac=cfv.scalfact2(ex1,ey1,es1[:,0],0.2) +cfv.secforce2(ex1,ey1,es1[:,0],plotpar,sfac) +cfv.secforce2(ex2,ey2,es2[:,0],plotpar,sfac) +cfv.secforce2(ex3,ey3,es3[:,0],plotpar,sfac) +cfv.plt.axis([-1.5, 7.5, -0.5, 5.5]) +cfv.title("Normal force") + +#----- Draw shear force diagram --------------------------------- +cfv.figure(3,fig_size=(6,4)) +plotpar=[2,1] +sfac=cfv.scalfact2(ex3,ey3,es3[:,1],0.2) +cfv.secforce2(ex1,ey1,es1[:,1],plotpar,sfac) +cfv.secforce2(ex2,ey2,es2[:,1],plotpar,sfac) +cfv.secforce2(ex3,ey3,es3[:,1],plotpar,sfac) +cfv.plt.axis([-1.5, 7.5, -0.5, 5.5]) +cfv.title("Shear force") + +#----- Draw moment diagram -------------------------------------- +cfv.figure(4,fig_size=(6,4)) +plotpar=[2,1] +sfac=cfv.scalfact2(ex3,ey3,es3[:,2],0.2) +cfv.secforce2(ex1,ey1,es1[:,2],plotpar,sfac) +cfv.secforce2(ex2,ey2,es2[:,2],plotpar,sfac) +cfv.secforce2(ex3,ey3,es3[:,2],plotpar,sfac) +cfv.plt.axis([-1.5, 7.5, -0.5, 5.5]) +cfv.title("Bending moment") +cfv.showAndWait() +#------------------------ end ----------------------------------- \ No newline at end of file From 164b3af6a5100300e53629d778c0e07e8be3055a Mon Sep 17 00:00:00 2001 From: Adam Sciegaj Date: Sun, 18 Dec 2022 17:55:29 +0100 Subject: [PATCH 5/7] Added bar2gs function and examples exn_bar2g, exn_beam2, exn_beam2_b, enx_bar2m --- calfem/core.py | 61 +++++++++++++++++ examples/exn_bar2g.py | 67 ++++++++++++++++++ examples/exn_bar2m.py | 98 ++++++++++++++++++++++++++ examples/exn_beam2.py | 148 ++++++++++++++++++++++++++++++++++++++++ examples/exn_beam2_b.py | 112 ++++++++++++++++++++++++++++++ 5 files changed, 486 insertions(+) create mode 100644 examples/exn_bar2g.py create mode 100644 examples/exn_bar2m.py create mode 100644 examples/exn_beam2.py create mode 100644 examples/exn_beam2_b.py diff --git a/calfem/core.py b/calfem/core.py index ffa56c7..ac16b25 100644 --- a/calfem/core.py +++ b/calfem/core.py @@ -226,6 +226,67 @@ def bar2s(ex, ey, ep, ed): return N.item() +def bar2gs(ex, ey, ep, ed): + """ + Calculate section forces in a two dimensional geometric + nonlinear bar element (bar2g). + Parameters: + ex = [x1 x2] element node coordinates + ey = [y1 y2] + + ep = [E A] element properties; + E: Young's modulus + A: cross section area + + ed = [u1 ... u4] element displacement vector + + Returns: + es = [N1; + N2 ] section forces, local directions + + QX: axial force + + edi = [ u1 ; element displacements, local directions, + u2 ; in n points along the bar, dim(es)= n x 1 + ...] + + eci = [ x1 ; local x-coordinates of the evaluation + x2 ; points, (x1=0 and xn=L) + ...] + """ + EA = ep[0]*ep[1] + ne = 2 + + dx = ex[1] - ex[0] + dy = ey[1] - ey[0] + L = np.sqrt(dx**2 + dy**2) + + n = [dx/L, dy/L, -dy/L, dx/L] + G = np.array([ + [n[0], n[1], 0., 0.], + [n[2], n[3], 0., 0.], + [0., 0., n[0], n[1]], + [0., 0., n[2], n[3]] + ]) + + edl = G@ed.reshape(-1,1) + a1 = np.array([edl[0], edl[2]]).reshape(-1,1) + + C1 = np.array([[1., 0.,], + [-1/L, 1/L]]) + C1a = C1@a1 + + x = np.linspace(0,L,ne).reshape(-1,1) + zero = np.zeros(x.shape) + one = np.ones(x.shape) + + u = np.concatenate((one, x),axis=1)@C1a + du = np.concatenate((zero, one),axis=1)@C1a + + N = EA*du + return N, N[0].item(), u, x + + def bar3e(ex, ey, ez, ep): """ Compute element stiffness matrix for three dimensional bar element. diff --git a/examples/exn_bar2g.py b/examples/exn_bar2g.py new file mode 100644 index 0000000..d0ab094 --- /dev/null +++ b/examples/exn_bar2g.py @@ -0,0 +1,67 @@ +# example exn_bar2g +#-------------------------------------------------------------------------- +# PURPOSE +# Analysis of a plane truss using second order theory. +#-------------------------------------------------------------------------- + +import numpy as np +import calfem.core as cfc + +# ----- Topology ----- + +edof = np.array([[1, 2, 5, 6], + [3, 4, 5, 6]]) + +# ----- Element properties and global coordinates ------------------------- + +E=10e9 +A1=4e-2 +A2=1e-2 +ep1=np.array([E, A1]) +ep2=np.array([E, A2]) + +ex1=np.array([0., 1.6]) +ey1=np.array([0., 0.]) +ex2=np.array([0., 1.6]) +ey2=np.array([1.2, 0]) + +# ----- Initial values for the iteration ---------------------------------- + +eps=1e-6 # Error norm +QX1=0.01 +QX2=0; # Initial axial forces +QX01=1 # Axial force of the initial former iteration +n=0 # Iteration counter + +# ----- Iteration procedure ----------------------------------------------- +while abs((QX1-QX01)/QX01) > eps: + n += 1 + + K = np.zeros((6,6)) + f = np.zeros((6,1)) + f[4] = -10e6 + f[5] = -0.2e6 + + Ke1 = cfc.bar2g(ex1,ey1,ep1,QX1) + Ke2 = cfc.bar2g(ex2,ey2,ep2,QX2) + K = cfc.assem(edof[0,:],K,Ke1) + K=cfc.assem(edof[1,:],K,Ke2) + bc = np.array([1,2,3,4]) + a, r = cfc.solveq(K,f,bc) + + Ed = cfc.extract_ed(edof,a) + + QX01 = QX1 + es1, QX1, ed1, _ = cfc.bar2gs(ex1,ey1,ep1,Ed[0,:]) + es2, QX2, ed2, _ = cfc.bar2gs(ex2,ey2,ep2,Ed[1,:]) + + if n>20: + print("The solution does not converge") + break + +# ----- Results ----------------------------------------------- +print("Displacements:") +print(a) +print("Normal forces:") +print(QX1) +print(QX2) \ No newline at end of file diff --git a/examples/exn_bar2m.py b/examples/exn_bar2m.py new file mode 100644 index 0000000..9e72afc --- /dev/null +++ b/examples/exn_bar2m.py @@ -0,0 +1,98 @@ +# example exn_bar2m +#-------------------------------------------------------------------------- +# PURPOSE +# Analysis of a plane truss considering material nonlinearity. +#-------------------------------------------------------------------------- + +import numpy as np +import calfem.core as cfc +import calfem.vis_mpl as cfv + +edof = np.array([ + [1, 2, 5, 6], + [5, 6, 7, 8], + [3, 4, 5, 6] +]) + +bc = np.array([1, 2, 3, 4, 7, 8]) + +ex = np.array([ + [0., 1.6], + [1.6, 1.6], + [0., 1.6] +]) + +ey = np.array([ + [0., 0.], + [0., 1.2], + [1.2, 0.] +]) + +E = np.array([200e9, 200e9, 200e9]) +A = np.array([6.0e-4, 3.0e-4, 10.0e-4]) +SY = 400e6 +Ns = (SY*A).reshape(-1,1) + +dp = 4e3 + +incr = 100 + +a = np.zeros((8,1)) +r = np.zeros((8,1)) +es = np.zeros((3,1)) + +plbar = 0 +pl = np.array([0., 0.]) + +# Forward Euler increamental solution + +for i in range(incr): + K = np.zeros((8,8)) + df = np.zeros((8,1)) + df[5] = -dp + + # Create and assemble element tangent stiffness matrix + for j in range(3): + ep = np.array([E[j], A[j]]) + Ke = cfc.bar2e(ex[j,:],ey[j,:],ep) + K = cfc.assem(edof[j,:], K, Ke) + + # Stop iteration if determinant det(Kr) <= 0 + fdof = np.setdiff1d(np.arange(1,9),bc) - 1 + Kr = K[np.ix_(fdof,fdof)] + if np.linalg.det(Kr) <= 0: + print("Determinant zero after increment ", i) + break + + # Solve for the displacement increment and determine total displacements + da, dr = cfc.solveq(K,df,bc) + a += da + r += dr + + # Determine normal forces in elements + ded = cfc.extract_ed(edof,da) + des = np.zeros((3,1)) + for j in range(3): + ep = np.array([E[j], A[j]]) + desj = cfc.bar2s(ex[j,:],ey[j,:],ep,ded[j,:]) + des[j,0] = desj + es += des + for j in range(3): + if abs(es[j,0]) >= Ns[j]: + E[j] = 0 + + # Determine if the stress in a bar has reached the yield stress + newplbar = np.sum(abs(es)>Ns) + if newplbar > plbar: + plbar = newplbar + print(plbar, "plastic elements for increment ", i+1, " at load = ", (i+1)*dp) + + # Save variables for curve plotting + pl = np.vstack((pl, np.array([-a[5,0], (i+1)*dp]))) + +# Plot force-displacement relation +cfv.figure(1, fig_size=(7,4)) +cfv.plt.plot(pl[:,0],pl[:,1]) +cfv.plt.xlabel("Displacement") +cfv.plt.ylabel("Force") +cfv.showAndWait() \ No newline at end of file diff --git a/examples/exn_beam2.py b/examples/exn_beam2.py new file mode 100644 index 0000000..4a8ad3e --- /dev/null +++ b/examples/exn_beam2.py @@ -0,0 +1,148 @@ +# example exn_beam2g +# ---------------------------------------------------------------- +# PURPOSE +# Geometrically nonlinear analysis of a plane frame. +# ---------------------------------------------------------------- + +import numpy as np +import calfem.core as cfc +import calfem.vis_mpl as cfv + + +#----- Topology ------------------------------------------------- + +edof = np.array([ + [4, 5, 6, 1, 2, 3], + [7, 8, 9, 10, 11, 12], + [4, 5, 6, 7, 8, 9] +]) + +#----- Element stiffness and element load matrices ------------- + +E=200e9 +A1=2e-3 +A2=6e-3 +I1=1.6e-5 +I2=5.4e-5 + +ep1 = np.array([E, A1, I1]) +ep3 = np.array([E, A2, I2]) + +ex1 = np.array([0.,0.]) +ey1 = np.array([4.,0.]) +ex2 = np.array([6.,6.]) +ey2 = np.array([4.,0.]) +ex3 = np.array([0.,6.]) +ey3 = np.array([4.,4.]) + +eq1 = np.array([0.]) +eq2 = np.array([0.]) +eq3 = np.array([-50e3]) + +# ----- Initial axial forces ---------------------------------------------- + +QX1 = 1e-4 +QX2 = 0 +QX3 = 0 +QX01 = 1 + +# ----- Iteration for convergence ----------------------------------------- + +eps = 1e-6 +n=0 + +while abs((QX1-QX01)/QX01) > eps: + n += 1 + K = np.zeros([12,12]) + f = np.zeros([12,1]) + f[3,0] = 10e3 + + Ke1 = cfc.beam2g(ex1,ey1,ep1,QX1) + Ke2 = cfc.beam2g(ex2,ey2,ep1,QX2) + Ke3, fe3 = cfc.beam2g(ex3,ey3,ep3,QX3,eq3) + + K = cfc.assem(edof[0,:],K,Ke1) + K = cfc.assem(edof[1,:],K,Ke2) + K, f = cfc.assem(edof[2,:],K,Ke3,f,fe3) + + bc = np.array([1, 2, 3, 10, 11]) + a, r = cfc.solveq(K, f, bc) + + ed = cfc.extract_ed(edof,a) + + QX01=QX1 + es1 = cfc.beam2gs(ex1,ey1,ep1,ed[0,:],QX1,eq1) + es2 = cfc.beam2gs(ex2,ey2,ep1,ed[1,:],QX2,eq2) + es3 = cfc.beam2gs(ex3,ey3,ep3,ed[2,:],QX3,eq3) + QX1 = es1[0,0] + QX2 = es2[0,0] + QX3 = es3[0,0] + + if n==1: + ed0 = ed + + if n>20: + print("The solution does not converge") + break + + +# #----- Section forces --------------------------------------- + +eq1 = np.array([0.,eq1.item()]) +eq2 = np.array([0.,eq2.item()]) +eq3 = np.array([0.,eq3.item()]) +es1, edi1, eci1 = cfc.beam2s(ex1,ey1,ep1,ed[0,:],eq1,21) +es2, edi2, eci2 = cfc.beam2s(ex2,ey2,ep1,ed[1,:],eq2,21) +es3, edi3, eci3 = cfc.beam2s(ex3,ey3,ep3,ed[2,:],eq3,21) + + +#----- Draw deformed frame --------------------------------------- +cfv.figure(1,fig_size=(6,4)) +plotpar=[3,1,0] +cfv.eldraw2(ex1,ey1,plotpar) +cfv.eldraw2(ex2,ey2,plotpar) +cfv.eldraw2(ex3,ey3,plotpar) +sfac=cfv.scalfact2(ex3,ey3,edi3,0.1) +plotpar=[1,2,1] +cfv.eldisp2(ex1,ey1,ed[0,:],plotpar,sfac) +cfv.eldisp2(ex2,ey2,ed[1,:],plotpar,sfac) +cfv.eldisp2(ex3,ey3,ed[2,:],plotpar,sfac) +plotpar=[2,4,2] +cfv.eldisp2(ex1,ey1,ed0[0,:],plotpar,sfac) +cfv.eldisp2(ex2,ey2,ed0[1,:],plotpar,sfac) +cfv.eldisp2(ex3,ey3,ed0[2,:],plotpar,sfac) +cfv.plt.axis([-1.5, 7.5, -0.5, 5.5]) +cfv.title("Displacements") + + +#----- Draw normal force diagram -------------------------------- +cfv.figure(2,fig_size=(6,4)) +plotpar=[2,1] +sfac=cfv.scalfact2(ex1,ey1,es1[:,0],0.2) +cfv.secforce2(ex1,ey1,es1[:,0],plotpar,sfac) +cfv.secforce2(ex2,ey2,es2[:,0],plotpar,sfac) +cfv.secforce2(ex3,ey3,es3[:,0],plotpar,sfac) +cfv.plt.axis([-1.5, 7.5, -0.5, 5.5]) +cfv.title("Normal force") + +#----- Draw shear force diagram --------------------------------- +cfv.figure(3,fig_size=(6,4)) +plotpar=[2,1] +sfac=cfv.scalfact2(ex3,ey3,es3[:,1],0.2) +cfv.secforce2(ex1,ey1,es1[:,1],plotpar,sfac) +cfv.secforce2(ex2,ey2,es2[:,1],plotpar,sfac) +cfv.secforce2(ex3,ey3,es3[:,1],plotpar,sfac) +cfv.plt.axis([-1.5, 7.5, -0.5, 5.5]) +cfv.title("Shear force") + +#----- Draw moment diagram -------------------------------------- +cfv.figure(4,fig_size=(6,4)) +plotpar=[2,1] +sfac=cfv.scalfact2(ex3,ey3,es3[:,2],0.2) +cfv.secforce2(ex1,ey1,es1[:,2],plotpar,sfac) +cfv.secforce2(ex2,ey2,es2[:,2],plotpar,sfac) +cfv.secforce2(ex3,ey3,es3[:,2],plotpar,sfac) +cfv.plt.axis([-1.5, 7.5, -0.5, 5.5]) +cfv.title("Bending moment") +cfv.showAndWait() +#------------------------ end ----------------------------------- \ No newline at end of file diff --git a/examples/exn_beam2_b.py b/examples/exn_beam2_b.py new file mode 100644 index 0000000..471b301 --- /dev/null +++ b/examples/exn_beam2_b.py @@ -0,0 +1,112 @@ +# example exn_beam2g_b +# ---------------------------------------------------------------- +# PURPOSE +# Buckling analysis of a plane frame. +# ---------------------------------------------------------------- + +import numpy as np +import calfem.core as cfc +import calfem.vis_mpl as cfv + + +#----- Topology ------------------------------------------------- + +edof = np.array([ + [4, 5, 6, 1, 2, 3], + [7, 8, 9, 10, 11, 12], + [4, 5, 6, 7, 8, 9] +]) + +#----- Element stiffness and element load matrices ------------- + +E=200e9 +A1=2e-3 +A2=6e-3 +I1=1.6e-5 +I2=5.4e-5 + +ep1 = np.array([E, A1, I1]) +ep3 = np.array([E, A2, I2]) + +ex1 = np.array([0.,0.]) +ey1 = np.array([4.,0.]) +ex2 = np.array([6.,6.]) +ey2 = np.array([4.,0.]) +ex3 = np.array([0.,6.]) +ey3 = np.array([4.,4.]) + +eq1 = np.array([0.]) +eq2 = np.array([0.]) +eq3 = np.array([-50e3]) + +# ----- Initial axial forces ---------------------------------------------- + +QX1 = 1e-4 +QX2 = 0 +QX3 = 0 +QX01 = 1 + +# ----- Iteration for convergence ----------------------------------------- + +eps = 1e-6 +n=0 + +while abs((QX1-QX01)/QX01) > eps: + n += 1 + K = np.zeros([12,12]) + f = np.zeros([12,1]) + f[3,0] = 10e3 + + Ke1 = cfc.beam2g(ex1,ey1,ep1,QX1) + Ke2 = cfc.beam2g(ex2,ey2,ep1,QX2) + Ke3, fe3 = cfc.beam2g(ex3,ey3,ep3,QX3,eq3) + + K = cfc.assem(edof[0,:],K,Ke1) + K = cfc.assem(edof[1,:],K,Ke2) + K, f = cfc.assem(edof[2,:],K,Ke3,f,fe3) + if n==1: + K0 = K + + bc = np.array([1, 2, 3, 10, 11]) + a, r = cfc.solveq(K, f, bc) + + ed = cfc.extract_ed(edof,a) + + QX01=QX1 + es1 = cfc.beam2gs(ex1,ey1,ep1,ed[0,:],QX1,eq1) + es2 = cfc.beam2gs(ex2,ey2,ep1,ed[1,:],QX2,eq2) + es3 = cfc.beam2gs(ex3,ey3,ep3,ed[2,:],QX3,eq3) + QX1 = es1[0,0] + QX2 = es2[0,0] + QX3 = es3[0,0] + + if n>20: + print("The solution does not converge") + break + + +# ----- Buckling analysis ------------------------------------------------- + +lam, phi = cfc.eigen(K,K0,bc) +one = np.ones(lam.shape) +alpha = np.divide(one,one-lam) +print(alpha[0]) + +# ----- Draw shape at instability ----------------------------------------- + +Ed = cfc.extract_ed(edof,-phi[:,0]) +cfv.figure(1,fig_size=(6,4)) +plotpar=[3,1,0] +cfv.eldraw2(ex1,ey1,plotpar) +cfv.eldraw2(ex2,ey2,plotpar) +cfv.eldraw2(ex3,ey3,plotpar) +sfac=cfv.scalfact2(ex3,ey3,Ed[2,:],0.1) +plotpar=[1,2,1] +cfv.eldisp2(ex1,ey1,Ed[0,:],plotpar,sfac) +cfv.eldisp2(ex2,ey2,Ed[1,:],plotpar,sfac) +cfv.eldisp2(ex3,ey3,Ed[2,:],plotpar,sfac) +cfv.plt.axis([-1.5, 7.5, -0.5, 5.5]) +cfv.title("Shape at instability") +cfv.showAndWait() + +#------------------------ end ----------------------------------- \ No newline at end of file From 01a2350431a2418cd9b810cf6e61c7979b1af3f3 Mon Sep 17 00:00:00 2001 From: Adam Sciegaj Date: Mon, 19 Dec 2022 11:10:26 +0100 Subject: [PATCH 6/7] Fixed a small bug in step2 --- calfem/core.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/calfem/core.py b/calfem/core.py index ac16b25..d821851 100644 --- a/calfem/core.py +++ b/calfem/core.py @@ -4073,7 +4073,7 @@ def step2(K,C,M,f,a0,da0,bc,ip,times,dofs): dim(dofhist['d2ahist']) = ndof x (nstep + 1) """ ndof, _ = K.shape - if not C: + if not np.array(C).any(): C = np.zeros((ndof,ndof)) dt, tottime, alpha, delta = ip b1 = dt*dt*0.5*(1-2*alpha) @@ -4108,7 +4108,7 @@ def step2(K,C,M,f,a0,da0,bc,ip,times,dofs): modelhist = {} sa=0 - if not times.any(): + if not np.array(times).any(): ntimes=0 sa=1 modelhist['a'] = np.zeros((ndof,nstep+1)) @@ -4123,7 +4123,7 @@ def step2(K,C,M,f,a0,da0,bc,ip,times,dofs): modelhist['d2a'] = np.zeros((ndof,ntimes)) dofhist = {} - if dofs.all(): + if np.array(dofs).all(): ndofs = len(dofs) if ndofs: dofhist['a'] = np.zeros((ndofs,nstep+1)) From 1615ee3010ae21fed20944448ef4f99fd8897189 Mon Sep 17 00:00:00 2001 From: Adam Sciegaj Date: Mon, 19 Dec 2022 12:46:22 +0100 Subject: [PATCH 7/7] Added step1 function to core --- calfem/core.py | 186 ++++++++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 183 insertions(+), 3 deletions(-) diff --git a/calfem/core.py b/calfem/core.py index d821851..6d937a8 100644 --- a/calfem/core.py +++ b/calfem/core.py @@ -4018,6 +4018,185 @@ def gfunc(G,dt): return ti, g1 +def step1(K,C,f,a0,bc,ip,times,dofs): + """ + Algorithm for dynamic solution of first-order + FE equations considering boundary conditions. + + Parameters: + + K conductivity matrix, dim(K) = ndof x ndof + C capacity matrix, dim(C) = ndof x ndof + f load vector, dim(f) = ndof x (nstep + 1), + If dim(f) = ndof x 1, the values are kept constant + during time integration + a0 initial vector a(0), dim(a0) = ndof x 1 + bc boundary condition matrix, dim(bc) = nbc x (nstep + 2) + where nbc = number of prescribed degrees of freedom (either constant or time-dependent) + The first column contains the numbers of the prescribed degrees of freedom + and the subsequent columns contain the time history. + If dim(bc) = nbc x 2, the values from the second column are kept constant + during time integration + ip array [dt, tottime, alpha], where + dt is the size of the time increment, + tottime is the total time, + alpha is time integration constant. + Frequently used values of alpha are: + alpha=0: forward difference; forward Euler, + alpha=1/2: trapezoidal rule; Crank-Nicholson + alpha=1: backward difference; backward Euler + times array [t(i) ...] of times at which output should be written to a and da + dofs array [dof(i) ...] of degree of freedom numbers for which history output + should be written to ahist and dahist + + Returns: + + modelhist dictionary containing solution history for the whole model at following keys: + modelhist['a'] constains values of a at all timesteps, + alternatively at times specified in 'times' + dim(modelhist['a']) = ndof x (nstep + 1) or ndof x ntimes + modelhist['da'] constains values of da at all timesteps, + alternatively at times specified in 'times' + dim(modelhist['da']) = ndof x (nstep + 1) or ndof x ntimes + dofhist dictionary containing solution history for the degrees of freedom selected in 'dofs': + dofhist['a'] constains time history of a at the dofs specified in 'dofs' + dim(dofhist['ahist']) = ndof x (nstep + 1) + dofhist['da'] constains time history of daat the dofs specified in 'dofs' + dim(dofhist['dahist']) = ndof x (nstep + 1) + """ + ndof, _ = K.shape + dt, tottime, alpha = ip + a1 = (1-alpha)*dt + a2 = alpha*dt + + nstep = 1 + if np.array(f).any(): + _, ncf = f.shape + if ncf>1: + nstep = ncf-1 + + if np.array(bc).any(): + _, ncb = bc.shape + if ncb>2: + nstep = ncb-2 + bound = 1 + if not np.array(bc).any(): + bound = 0 + + ns = int(tottime/dt) + if (ns < nstep or nstep==1): + nstep=ns + + tf = np.zeros((ndof,nstep+1)) + if np.array(f).any(): + if ncf==1: + tf = f[:,0].reshape(-1,1)@np.ones((1,nstep+1)) + if ncf>1: + tf = np.copy(f) + + modelhist = {} + sa=0 + if not np.array(times).any(): + ntimes=0 + sa=1 + modelhist['a'] = np.zeros((ndof,nstep+1)) + modelhist['da'] = np.zeros((ndof,nstep+1)) + else: + ntimes = len(times) + if ntimes: + sa=2 + modelhist['a'] = np.zeros((ndof,ntimes)) + modelhist['da'] = np.zeros((ndof,ntimes)) + + dofhist = {} + if np.array(dofs).all(): + ndofs = len(dofs) + if ndofs: + dofhist['a'] = np.zeros((ndofs,nstep+1)) + dofhist['da'] = np.zeros((ndofs,nstep+1)) + else: + ndofs=0 + + itime = 0 + + # Calculate initial second time derivative d2a0 + da0 = np.linalg.solve(C,tf[:,0].reshape(-1,1) - K@a0) + # Save initial values + if sa==1: + modelhist['a'][:,0] = a0.ravel() + modelhist['da'][:,0] = da0.ravel() + elif sa==2: + if times[itime]==0: + modelhist['a'][:,itime] = a0.ravel() + modelhist['da'][:,itime] = da0.ravel() + itime += 1 + + if ndofs: + dofhist['a'][:,0] = a0[np.ix_(dofs-1)].ravel() + dofhist['da'][:,0] = da0[np.ix_(dofs-1)].ravel() + + # Reduce matrices due to bcs + tempa = np.zeros((ndof,1)) + tempda = np.zeros((ndof,1)) + fdof=np.arange(1,ndof+1).astype(int) + if bound: + nrb, ncb = bc.shape + if ncb==2: + pa = bc[:,1].reshape(-1,1)@np.ones((1,nstep+1)) + pda = np.zeros((nrb,nstep+1)) + elif ncb>2: + pa = np.copy(bc[:,1:]) + pda1 = (pa[:,1]-pa[:,0])/dt + pdarest = (pa[:,1:] - pa[:,0:-1])/dt + pda = np.hstack((pda1.reshape(-1,1),pdarest)) + pdof = np.copy(bc[:,0]).astype(int) + fdof = np.setdiff1d(fdof,pdof).astype(int) - 1 + pdof -= 1 #adjusting for indexing starting from 0 + Keff = C[np.ix_(fdof,fdof)] + a2*K[np.ix_(fdof,fdof)] + else: + fdof -= 1 #adjusting for indexing starting from 0 + Keff = C + a2*K + + L, U = lu(Keff,permute_l=True) + anew = a0[np.ix_(fdof)] + danew = da0[np.ix_(fdof)] + + # Iterate over time steps + for j in range(1,nstep+1): + time = dt*j + aold = np.copy(anew) + daold = np.copy(danew) + apred = aold + a1*daold + if not bound: + reff = tf[:,j].reshape(-1,1) - K@apred + else: + pdeff = C[np.ix_(fdof,pdof)]@pda[:,j].reshape(-1,1) + K[np.ix_(fdof,pdof)]@pa[:,j].reshape(-1,1) + reff = tf[np.ix_(fdof),j].reshape(-1,1) - K[np.ix_(fdof,fdof)]@apred - pdeff + y = np.linalg.solve(L,reff) + danew = np.linalg.solve(U,y) + anew = apred + a2*danew + # Save to modelhist and dofhist + if bound: + tempa[np.ix_(pdof)] = pa[:,j].reshape(-1,1) + tempda[np.ix_(pdof)] = pda[:,j].reshape(-1,1) + tempa[np.ix_(fdof)] = anew + tempda[np.ix_(fdof)] = danew + if sa==1: + modelhist['a'][:,j] = tempa.ravel() + modelhist['da'][:,j] = tempda.ravel() + elif sa==2: + if ntimes and itime < ntimes: + if time >= times[itime]: + modelhist['a'][:,itime] = tempa.ravel() + modelhist['da'][:,itime] = tempda.ravel() + itime += 1 + if ndofs: + dofhist['a'][:,j] = tempa[np.ix_(dofs-1)].ravel() + dofhist['da'][:,j] = tempda[np.ix_(dofs-1)].ravel() + + return modelhist, dofhist + + def step2(K,C,M,f,a0,da0,bc,ip,times,dofs): """ Algorithm for dynamic solution of second-order @@ -4040,11 +4219,11 @@ def step2(K,C,M,f,a0,da0,bc,ip,times,dofs): and the subsequent columns contain the time history. If dim(bc) = nbc x 2, the values from the second column are kept constant during time integration - ip list [dt, tottime, alpha, delta], where + ip array [dt, tottime, alpha, delta], where dt is the size of the time increment, tottime is the total time, alpha and delta are time integration constants for the Newmark family of methods. - Frequently used values of alpha and delta: + Frequently used values of alpha and delta are: alpha=1/4, delta=1/2: average acceleration (trapezoidal) rule, alpha=1/6, delta=1/2: linear acceleration alpha=0, delta=1/2: central difference @@ -4064,7 +4243,7 @@ def step2(K,C,M,f,a0,da0,bc,ip,times,dofs): modelhist['d2a'] constains acceleration values at all timesteps, alternatively at times specified in 'times' dim(modelhist['d2a']) = ndof x (nstep + 1) or ndof x ntimes - dofhist dictionary containing solution history for the selected 'dofs': + dofhist dictionary containing solution history for the degrees of freedom selected in 'dofs': dofhist['a'] constains displacement time history at the dofs specified in 'dofs' dim(dofhist['ahist']) = ndof x (nstep + 1) dofhist['da'] constains velocity time history at the dofs specified in 'dofs' @@ -4223,6 +4402,7 @@ def step2(K,C,M,f,a0,da0,bc,ip,times,dofs): return modelhist, dofhist + def extract_eldisp(edof, a): """ Extract element displacements from the global displacement