Added new commands and and updated links

calfem/core.py

@@ -6,6 +6,7 @@
 """
 
 from scipy.sparse.linalg import dsolve
+from scipy.linalg import eig
 import numpy as np
 
 import logging as cflog
@@ -298,6 +299,113 @@ def bar3s(ex, ey, ez, ep, ed):
 
     return N.item()
 
+def beam1e(ex, ep, eq=None):
+    """
+    Compute the stiffness matrix for a one dimensional beam element.
+
+    :param list ex: element x coordinates [x1, x2]
+    :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)
+    """
+    L = ex[1]-ex[0]
+
+    E = ep[0]
+    I = ep[1]
+
+    qy = 0.
+    if not eq is None:
+        qy = eq
+
+    Ke = E*I/(L**3) * np.mat([
+        [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**2]
+    ])
+
+    fe = qy*np.mat([L/2, L**2/12, L/2, -L**2/12]).T
+
+    if eq is None:
+        return Ke
+    else:
+        return Ke, fe
+
+def beam1s(ex, ep, ed, eq=None, nep=None):
+    """
+    Compute section forces in one dimensional beam element (beam1e).
+
+    Parameters:
+
+        ex = [x1 x2]        element node coordinates
+
+        ep = [E I]          element properties,
+                            E:  Young's modulus
+                            I:  moment of inertia
+
+        ed = [u1 ... u4]    element displacements
+
+        eq = qy             distributed load, local directions
+
+        nep                 number of evaluation points ( default=2 )
+
+    Returns:
+
+        es = [ V1 M1        section forces, local directions, in
+               V2 M2        n points along the beam, dim(es)= n x 2
+               .........]
+
+        edi = [ v1          element displacements, local directions,
+                v2          in n points along the beam, dim(es)= n x 1
+                .......]
+
+            eci = [ x1      local x-coordinates of the evaluation
+                    x2      points, (x1=0 and xn=L)
+                    ...]
+
+    """
+    EI = ep[0]*ep[1]
+    L = ex[1]-ex[0]
+
+    qy = 0.
+
+    if not eq is None:
+        qy = eq
+
+    ne = 2
+
+    if nep != None:
+        ne = nep
+
+    Cinv = np.mat([
+        [1, 0,  0,  0],
+        [0, 1,  0,  0],
+        [-3/(L**2), -2/L,   3/(L**2),   -1/L],
+        [2/(L**3),  1/(L**2),   -2/(L**3),  1/(L**2)]
+    ])
+
+    Ca = (Cinv@ed).T
+
+    x = np.asmatrix(np.arange(0., L+L/(ne-1), L/(ne-1))).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)
+
+    M = EI*d2v
+    V = -EI*d3v
+    edi = v
+    eci = x
+    es = np.concatenate((V, M), 1)
+
+    return (es, edi, eci)
+
+
 
 def beam2e(ex, ey, ep, eq=None):
     """
@@ -3781,6 +3889,53 @@ def spsolveq(K, f, bcPrescr, bcVal=None):
     return (a_m, Q)
 
 
+def eigen(K,M,b=None):
+    """
+    Solve the generalized eigenvalue problem
+    |K-LM|X = 0, considering boundary conditions
+
+    Parameters:
+
+        K           global stiffness matrix, dim(K) = ndof x ndof
+        M           global mass matrix, dim(M) = ndof x ndof
+        b           boundary condition vector, dim(b) = nbc x 1
+
+    Returns:
+
+        L           eigenvalue vector, dim(L) = (ndof-nbc) x 1
+        X           eigenvectors, dim(X) = ndof x (ndof-nbc)
+    """
+    nd, _ = K.shape
+    if b is not None:
+        fdof = np.setdiff1d(np.arange(nd), b-1)
+        D, X1 = eig(K[np.ix_(fdof,fdof)], M[np.ix_(fdof,fdof)])
+        D = np.real(D)
+        nfdof, _ = X1.shape
+        for j in range(nfdof):
+            mnorm = np.sqrt(X1[:,j].T@M[np.ix_(fdof,fdof)]@X1[:,j])
+            X1[:,j] /= mnorm
+        s_order = np.argsort(D)
+        L = np.sort(D)
+        X2 = np.zeros(X1.shape)
+        for ind,j in enumerate(s_order):
+            X2[:,ind] = X1[:,j]
+        X = np.zeros((nd,nfdof))
+        X[fdof,:] = X2
+        return L, X
+    else:
+        D, X1 = eig(K, M)
+        D = np.real(D)
+        for j in range(nd):
+            mnorm = np.sqrt(X1[:,j].T@M@X1[:,j])
+            X1[:,j] /= mnorm
+        s_order = np.argsort(D)
+        L = np.sort(D)
+        X = np.copy(X1)
+        for ind,j in enumerate(s_order):
+            X[:,ind] = X1[:,j]
+        return L, X
+
+
 def extract_eldisp(edof, a):
     """
     Extract element displacements from the global displacement

calfem/vis_mpl.py

@@ -1580,3 +1580,128 @@ def eldisp2(ex, ey, ed, plotpar=[2, 1, 1], sfac=None):
 #     plot(x,y,s2)
 #     hold off
 # %--------------------------end--------------------------------
+
+def secforce2(ex, ey, es, plotpar=None, sfac=None, eci = None):
+    """
+    --------------------------------------------------------------------------
+     PURPOSE:
+      Draw section force diagram for a two dimensional bar or beam element.
+
+     INPUT:
+        ex = [ x1 x2 ]
+        ey = [ y1 y2 ]	element node coordinates.
+
+        es = [ S1;
+               S2;
+             ... ] 	vector containing the section force
+                    in Nbr evaluation points along the element.
+
+        plotpar=[linecolour, elementcolour]
+
+                linecolour=1 -> black     elementcolour=1 -> black
+                           2 -> blue                    2 -> blue
+                           3 -> magenta                 3 -> magenta
+                           4 -> red                     4 -> red
+
+        sfac = [scalar]	scale factor for section force diagrams.
+
+        eci = [  x1;
+                x2;
+               ... ]  local x-coordinates of the evaluation points (Nbr).
+           If not given, the evaluation points are assumed to be uniformly
+           distributed
+    """
+
+    if ex.shape == ey.shape:
+        if ex.ndim !=1:
+            nen = ex.shape[1]
+        else:
+            nen = ex.shape[0]
+
+            ex = ex.reshape(1, nen)
+            ey = ey.reshape(1, nen)
+    else:
+        raise ValueError("Check size of ex, ey dimensions.")
+
+    #bar2s returns a float - convert it to array
+    if isinstance(es, float):
+        es = np.array([es,es]).reshape(2,1)
+
+    Nbr = es.shape[0]
+    b = np.array([ex[0,1]-ex[0,0],ey[0,1]-ey[0,0]])
+    Length = np.sqrt(b@b)
+    n = (b/Length).reshape(1,len(b))
+
+    if plotpar is None:
+        mpl_diagram_color = (0, 0, 1)  # 'b'
+        mpl_elem_color = (0, 0, 0)  # 'k'
+    elif len(plotpar) != 2:
+        raise ValueError('Check size of "plotpar" input argument!')
+    else:
+        p1, p2 = plotpar
+        if p1 == 1:
+            mpl_diagram_color = (0, 0, 0) # 'k'
+        elif p1 == 2:
+            mpl_diagram_color = (0, 0, 1) # 'b'
+        elif p1 == 3:
+            mpl_diagram_color = (1, 0, 1) # 'm'
+        elif p1 == 4:
+            mpl_diagram_color = (1, 0, 0) # 'r'
+
+        if p2 == 1:
+            mpl_elem_color = (0, 0, 0) # 'k'
+        elif p2 == 2:
+            mpl_elem_color = (0, 0, 1) # 'b'
+        elif p2 == 3:
+            mpl_elem_color = (1, 0, 1) # 'm'
+        elif p2 == 4:
+            mpl_elem_color = (1, 0, 0) # 'r'
+
+
+    if sfac is None:
+        sfac=(0.2*Length)/np.max(np.abs(es))
+
+    if eci is None:
+        eci = np.linspace(0,Length,Nbr).reshape(Nbr,1)
+
+    if es.shape[0] != eci.shape[0]:
+        raise ValueError("Check size of 'es' or 'eci' input argument!")
+
+    es *= sfac
+
+    # From local x-coordinates to global coordinates of the element
+    A = np.zeros([Nbr,2])
+    A = [ex[0][0], ey[0][0]] + eci@n
+    Ax = np.zeros([Nbr-1,2])
+    Ax[:,0] = A[:-1,0]
+    Ax[:,1] = A[1:,0]
+    Ay = np.zeros([Nbr-1,2])
+    Ay[:,0] = A[:-1,1]
+    Ay[:,1] = A[1:,1]
+
+
+    Bx=np.copy(Ax)
+    By=np.copy(Ay)
+    for i in range(Nbr-1):
+        Ax[i,:] = [ Ax[i,0]+es[i]*n[0,1], Ax[i,1]+es[i+1]*n[0,1] ]
+        Ay[i,:] = [ Ay[i,0]-es[i]*n[0,0], Ay[i,1]-es[i+1]*n[0,0] ]
+
+    #plot diagram and its vertical stripes at the ends
+    plt.axis('equal')
+    draw_elements(Ax, Ay, color=mpl_diagram_color,
+                  line_style='solid', filled=False, closed=False)
+
+    draw_elements(np.array([ex[0][0], Ax[0,0]]).reshape(1,2), np.array([ey[0][0], Ay[0,0]]).reshape(1,2), color=mpl_diagram_color,
+                  line_style='solid', filled=False, closed=False)
+    draw_elements(np.array([ex[0][1], Ax[-1,-1]]).reshape(1,2), np.array([ey[0][1], Ay[-1,-1]]).reshape(1,2), color=mpl_diagram_color,
+                  line_style='solid', filled=False, closed=False)
+
+    #plot stripes in diagram
+    for i in range(0,Nbr-2):
+        np.array([Bx[i,1], Ax[i,1]]).reshape(1,2), np.array([By[i,1], Ay[i,1]]).reshape(1,2)
+        draw_elements(np.array([Bx[i,1], Ax[i,1]]).reshape(1,2), np.array([By[i,1], Ay[i,1]]).reshape(1,2), color=mpl_diagram_color,
+                  line_style='solid', filled=False, closed=False)
+
+    #plot elements
+    draw_elements(ex, ey, color=mpl_elem_color,
+                  line_style='solid', filled=False, closed=False)