From c8f2a1b00e56cd3ddf216379224abe2e9e901fb4 Mon Sep 17 00:00:00 2001 From: Brooks Smith Date: Sat, 7 Oct 2023 19:27:55 +1100 Subject: [PATCH] Improve typing and formatting --- pycufsm/analysis_c.pyx | 1122 +++++++++++++++++++++++++--------------- pycufsm/analysis_p.py | 117 ++++- pycufsm/cfsm.py | 45 +- 3 files changed, 827 insertions(+), 457 deletions(-) diff --git a/pycufsm/analysis_c.pyx b/pycufsm/analysis_c.pyx index 80afe5f..c55afca 100644 --- a/pycufsm/analysis_c.pyx +++ b/pycufsm/analysis_c.pyx @@ -45,7 +45,7 @@ cpdef int constr_bc_flag(np.ndarray nodes, np.ndarray constraints): Returns: bc_flag (int): 1 if there are user constraints or node fixities 0 if there is no user constraints and node fixities - + Z. Li, June 2010 """ # Check for boundary conditions on the nodes @@ -60,8 +60,7 @@ cpdef int constr_bc_flag(np.ndarray nodes, np.ndarray constraints): # Check for user defined constraints too if len(constraints) == 0: return 0 - else: - return 1 + return 1 cpdef np.ndarray elem_prop(np.ndarray nodes, np.ndarray elements): @@ -106,20 +105,25 @@ cpdef np.ndarray elem_prop(np.ndarray nodes, np.ndarray elements): cpdef k_kg_global( - np.ndarray nodes, np.ndarray elements, np.ndarray el_props, np.ndarray props, - double length, str b_c, np.ndarray m_a + np.ndarray nodes, + np.ndarray elements, + np.ndarray el_props, + np.ndarray props, + double length, + str b_c, + np.ndarray m_a, ): """Generates element stiffness matrix (k_global) in global coordinates Generate geometric stiffness matrix (kg_global) in global coordinates Args: - nodes (np.ndarray): - elements (np.ndarray): - el_props (np.ndarray): - props (np.ndarray): - length (np.ndarray): - b_c (str): - m_a (np.ndarray): + nodes (np.ndarray): + elements (np.ndarray): + el_props (np.ndarray): + props (np.ndarray): + length (np.ndarray): + b_c (str): + m_a (np.ndarray): Returns: k_global (np.ndarray): global stiffness matrix @@ -187,7 +191,7 @@ cpdef k_kg_global( ty_2=ty_2, b_strip=b_strip, b_c=b_c, - m_a=m_a + m_a=m_a, ) # Transform k_local and kg_local into global coordinates @@ -205,15 +209,25 @@ cpdef k_kg_global( node_i=node_i, node_j=node_j, n_nodes=n_nodes, - m_a=m_a + m_a=m_a, ) return k_global, kg_global cdef k_kg_local( - double stiff_x, double stiff_y, double nu_x, double nu_y, double bulk, double thick, - double length, double ty_1, double ty_2, double b_strip, str b_c, np.ndarray m_a + double stiff_x, + double stiff_y, + double nu_x, + double nu_y, + double bulk, + double thick, + double length, + double ty_1, + double ty_2, + double b_strip, + str b_c, + np.ndarray m_a, ): """Generate element stiffness matrix (k_local) in local coordinates Generate geometric stiffness matrix (kg_local) in local coordinates @@ -239,10 +253,10 @@ cdef k_kg_local( m_a (np.ndarray): longitudinal terms (or half-wave numbers) for this length Returns: - k_local (np.ndarray): local stiffness matrix, a total_m x total_m matrix of 8 by 8 + k_local (np.ndarray): local stiffness matrix, a total_m x total_m matrix of 8 by 8 submatrices. k_local=[k_mp]total_m x total_m block matrix each k_mp is the 8 x 8 submatrix in the DOF order [u1 v1 u2 v2 w1 theta1 w2 theta2]' - kg_local (np.ndarray): local geometric stiffness matrix, a total_m x total_m matrix of + kg_local (np.ndarray): local geometric stiffness matrix, a total_m x total_m matrix of 8 by 8 submatrices. kg_local=[kg_mp]total_m x total_m block matrix each kg_mp is the 8 x 8 submatrix in the DOF order [u1 v1 u2 v2 w1 theta1 @@ -251,11 +265,11 @@ cdef k_kg_local( modified by Z. Li, June 2010 klocal and kglocal merged by B Smith, May 2023 """ - cdef double e_1 = stiff_x / (1 - nu_x*nu_y) - cdef double e_2 = stiff_y / (1 - nu_x*nu_y) - cdef double d_x = stiff_x * thick**3 / (12 * (1 - nu_x*nu_y)) - cdef double d_y = stiff_y * thick**3 / (12 * (1 - nu_x*nu_y)) - cdef double d_1 = nu_x * stiff_y * thick**3 / (12 * (1 - nu_x*nu_y)) + cdef double e_1 = stiff_x / (1 - nu_x * nu_y) + cdef double e_2 = stiff_y / (1 - nu_x * nu_y) + cdef double d_x = stiff_x * thick**3 / (12 * (1 - nu_x * nu_y)) + cdef double d_y = stiff_y * thick**3 / (12 * (1 - nu_x * nu_y)) + cdef double d_1 = nu_x * stiff_y * thick**3 / (12 * (1 - nu_x * nu_y)) cdef double d_xy = bulk * thick**3 / 12 cdef int total_m = len(m_a) # Total number of longitudinal terms m @@ -280,7 +294,7 @@ cdef k_kg_local( [i_1, i_2, i_3, i_4, i_5] = bc_i1_5(b_c=b_c, m_i=m_a[i], m_j=m_a[j], length=length) - k_local[8 * i:8*i + 4, 8 * j:8*j + 4] = calc_km_mp( + k_local[8 * i : 8 * i + 4, 8 * j : 8 * j + 4] = calc_km_mp( e_1=e_1, e_2=e_2, c_1=c_1, @@ -293,9 +307,9 @@ cdef k_kg_local( i_2=i_2, i_3=i_3, i_4=i_4, - i_5=i_5 + i_5=i_5, ) - k_local[8*i + 4:8 * (i+1), 8*j + 4:8 * (j+1)] = calc_kf_mp( + k_local[8 * i + 4 : 8 * (i + 1), 8 * j + 4 : 8 * (j + 1)] = calc_kf_mp( d_x=d_x, d_y=d_y, d_1=d_1, @@ -305,10 +319,10 @@ cdef k_kg_local( i_2=i_2, i_3=i_3, i_4=i_4, - i_5=i_5 + i_5=i_5, ) - kg_local[8 * i:8*i + 4, 8 * j:8*j + 4] = calc_gm_mp( + kg_local[8 * i : 8 * i + 4, 8 * j : 8 * j + 4] = calc_gm_mp( u_i=u_i, u_j=u_j, b_strip=b_strip, @@ -316,9 +330,9 @@ cdef k_kg_local( ty_1=ty_1, ty_2=ty_2, i_4=i_4, - i_5=i_5 + i_5=i_5, ) - kg_local[8*i + 4:8 * (i+1), 8*j + 4:8 * (j+1)] = calc_gf_mp( + kg_local[8 * i + 4 : 8 * (i + 1), 8 * j + 4 : 8 * (j + 1)] = calc_gf_mp( ty_1=ty_1, ty_2=ty_2, b_strip=b_strip, i_5=i_5 ) @@ -326,8 +340,13 @@ cdef k_kg_local( cpdef np.ndarray kglobal_transv( - np.ndarray nodes, np.ndarray elements, np.ndarray el_props, np.ndarray props, - double length, str b_c, double m_i + np.ndarray nodes, + np.ndarray elements, + np.ndarray props, + double m_i, + double length, + str b_c, + np.ndarray el_props, ): """this routine creates the global stiffness matrix for planar displacements basically the same way as in the main program, however: @@ -394,7 +413,7 @@ cpdef np.ndarray kglobal_transv( length=length, b_strip=b_strip, b_c=b_c, - m_i=m_i + m_i=m_i, ) # Transform k_local and kg_local into global coordinates @@ -410,15 +429,23 @@ cpdef np.ndarray kglobal_transv( k_local=k_local, node_i=node_i, node_j=node_j, - n_nodes=n_nodes + n_nodes=n_nodes, ) return k_global_transv cdef np.ndarray klocal_transv( - double stiff_x, double stiff_y, double nu_x, double nu_y, double bulk, double thick, - double length, double b_strip, str b_c, double m_i + double stiff_x, + double stiff_y, + double nu_x, + double nu_y, + double bulk, + double thick, + double length, + double b_strip, + str b_c, + double m_i, ): """this routine creates the local stiffness matrix for bending terms basically the same way as in the main program, however: @@ -460,10 +487,10 @@ cdef np.ndarray klocal_transv( cdef double c_2 = u_p / length [i_1, _, _, _, _] = bc_i1_5(b_c=b_c, m_i=m_i, m_j=m_i, length=length) - cdef double i_2 = 0 - cdef double i_3 = 0 - cdef double i_4 = 0 - cdef double i_5 = 0 + cdef double i_2 = double(0.0) + cdef double i_3 = double(0.0) + cdef double i_4 = double(0.0) + cdef double i_5 = double(0.0) k_local[0:4, 0:4] = calc_km_mp( e_1=e_1, @@ -478,7 +505,7 @@ cdef np.ndarray klocal_transv( i_2=i_2, i_3=i_3, i_4=i_4, - i_5=i_5 + i_5=i_5, ) k_local[4:8, 4:8] = calc_kf_mp( d_x=d_x, @@ -490,14 +517,25 @@ cdef np.ndarray klocal_transv( i_2=i_2, i_3=i_3, i_4=i_4, - i_5=i_5 + i_5=i_5, ) return k_local cdef np.ndarray calc_km_mp( - double e_1, double e_2, double c_1, double c_2, double b_strip, double bulk, - double nu_x, double thick, double i_1, double i_2, double i_3, double i_4, double i_5 + double e_1, + double e_2, + double c_1, + double c_2, + double b_strip, + double bulk, + double nu_x, + double thick, + double i_1, + double i_2, + double i_3, + double i_4, + double i_5, ): """Calculate the membrane stiffness sub-matrix, used in the assembly of local stiffness matrices @@ -528,31 +566,40 @@ cdef np.ndarray calc_km_mp( cdef np.ndarray[np.double_t, ndim=2] km_mp = np.zeros((4, 4)) # assemble the matrix of Km_mp (membrane stiffness matrix) - km_mp[0, 0] = e_1*i_1/b_strip + bulk*b_strip*i_5/3 - km_mp[0, 1] = e_2 * nu_x * (-1 / 2 / c_2) * i_3 - bulk*i_5/2/c_2 - km_mp[0, 2] = -e_1 * i_1 / b_strip + bulk*b_strip*i_5/6 - km_mp[0, 3] = e_2 * nu_x * (-1 / 2 / c_2) * i_3 + bulk*i_5/2/c_2 - - km_mp[1, 0] = e_2 * nu_x * (-1 / 2 / c_1) * i_2 - bulk*i_5/2/c_1 - km_mp[1, 1] = e_2*b_strip*i_4/3/c_1/c_2 + bulk*i_5/b_strip/c_1/c_2 - km_mp[1, 2] = e_2 * nu_x * (1/2/c_1) * i_2 - bulk*i_5/2/c_1 - km_mp[1, 3] = e_2*b_strip*i_4/6/c_1/c_2 - bulk*i_5/b_strip/c_1/c_2 - - km_mp[2, 0] = -e_1 * i_1 / b_strip + bulk*b_strip*i_5/6 - km_mp[2, 1] = e_2 * nu_x * (1/2/c_2) * i_3 - bulk*i_5/2/c_2 - km_mp[2, 2] = e_1*i_1/b_strip + bulk*b_strip*i_5/3 - km_mp[2, 3] = e_2 * nu_x * (1/2/c_2) * i_3 + bulk*i_5/2/c_2 - - km_mp[3, 0] = e_2 * nu_x * (-1 / 2 / c_1) * i_2 + bulk*i_5/2/c_1 - km_mp[3, 1] = e_2*b_strip*i_4/6/c_1/c_2 - bulk*i_5/b_strip/c_1/c_2 - km_mp[3, 2] = e_2 * nu_x * (1/2/c_1) * i_2 + bulk*i_5/2/c_1 - km_mp[3, 3] = e_2*b_strip*i_4/3/c_1/c_2 + bulk*i_5/b_strip/c_1/c_2 + km_mp[0, 0] = e_1 * i_1 / b_strip + bulk * b_strip * i_5 / 3 + km_mp[0, 1] = e_2 * nu_x * (-1 / 2 / c_2) * i_3 - bulk * i_5 / 2 / c_2 + km_mp[0, 2] = -e_1 * i_1 / b_strip + bulk * b_strip * i_5 / 6 + km_mp[0, 3] = e_2 * nu_x * (-1 / 2 / c_2) * i_3 + bulk * i_5 / 2 / c_2 + + km_mp[1, 0] = e_2 * nu_x * (-1 / 2 / c_1) * i_2 - bulk * i_5 / 2 / c_1 + km_mp[1, 1] = e_2 * b_strip * i_4 / 3 / c_1 / c_2 + bulk * i_5 / b_strip / c_1 / c_2 + km_mp[1, 2] = e_2 * nu_x * (1 / 2 / c_1) * i_2 - bulk * i_5 / 2 / c_1 + km_mp[1, 3] = e_2 * b_strip * i_4 / 6 / c_1 / c_2 - bulk * i_5 / b_strip / c_1 / c_2 + + km_mp[2, 0] = -e_1 * i_1 / b_strip + bulk * b_strip * i_5 / 6 + km_mp[2, 1] = e_2 * nu_x * (1 / 2 / c_2) * i_3 - bulk * i_5 / 2 / c_2 + km_mp[2, 2] = e_1 * i_1 / b_strip + bulk * b_strip * i_5 / 3 + km_mp[2, 3] = e_2 * nu_x * (1 / 2 / c_2) * i_3 + bulk * i_5 / 2 / c_2 + + km_mp[3, 0] = e_2 * nu_x * (-1 / 2 / c_1) * i_2 + bulk * i_5 / 2 / c_1 + km_mp[3, 1] = e_2 * b_strip * i_4 / 6 / c_1 / c_2 - bulk * i_5 / b_strip / c_1 / c_2 + km_mp[3, 2] = e_2 * nu_x * (1 / 2 / c_1) * i_2 + bulk * i_5 / 2 / c_1 + km_mp[3, 3] = e_2 * b_strip * i_4 / 3 / c_1 / c_2 + bulk * i_5 / b_strip / c_1 / c_2 return km_mp * thick + cdef np.ndarray calc_kf_mp( - double d_x, double d_y, double d_1, double d_xy, double b_strip, double i_1, - double i_2, double i_3, double i_4, double i_5 + double d_x, + double d_y, + double d_1, + double d_xy, + double b_strip, + double i_1, + double i_2, + double i_3, + double i_4, + double i_5, ): """Calculate the flexural stiffness sub-matrix, used in the assembly of local stiffness matrices @@ -580,45 +627,159 @@ cdef np.ndarray calc_kf_mp( cdef np.ndarray[np.double_t, ndim=2] kf_mp = np.zeros((4, 4)) # assemble the matrix of Kf_mp (flexural stiffness matrix) - kf_mp[0, 0] = (5040*d_x*i_1 - 504*b_strip**2*d_1*i_2 - 504*b_strip**2*d_1*i_3 \ - + 156*b_strip**4*d_y*i_4 + 2016*b_strip**2*d_xy*i_5)/420/b_strip**3 - kf_mp[0, 1] = (2520*b_strip*d_x*i_1 - 462*b_strip**3*d_1*i_2 - 42*b_strip**3*d_1*i_3 \ - + 22*b_strip**5*d_y*i_4 + 168*b_strip**3*d_xy*i_5)/420/b_strip**3 - kf_mp[0, 2] = (-5040*d_x*i_1 + 504*b_strip**2*d_1*i_2 + 504*b_strip**2*d_1*i_3 \ - + 54*b_strip**4*d_y*i_4 - 2016*b_strip**2*d_xy*i_5)/420/b_strip**3 - kf_mp[0, 3] = (2520*b_strip*d_x*i_1 - 42*b_strip**3*d_1*i_2 - 42*b_strip**3*d_1*i_3 \ - - 13*b_strip**5*d_y*i_4 + 168*b_strip**3*d_xy*i_5)/420/b_strip**3 - - kf_mp[1, 0] = (2520*b_strip*d_x*i_1 - 462*b_strip**3*d_1*i_3 - 42*b_strip**3*d_1*i_2 \ - + 22*b_strip**5*d_y*i_4 + 168*b_strip**3*d_xy*i_5)/420/b_strip**3 - kf_mp[1, 1] = (1680*b_strip**2*d_x*i_1 - 56*b_strip**4*d_1*i_2 - 56*b_strip**4*d_1*i_3 \ - + 4*b_strip**6*d_y*i_4 + 224*b_strip**4*d_xy*i_5)/420/b_strip**3 - kf_mp[1, 2] = (-2520*b_strip*d_x*i_1 + 42*b_strip**3*d_1*i_2 + 42*b_strip**3*d_1*i_3 \ - + 13*b_strip**5*d_y*i_4 - 168*b_strip**3*d_xy*i_5)/420/b_strip**3 - kf_mp[1, 3] = (840*b_strip**2*d_x*i_1 + 14*b_strip**4*d_1*i_2 + 14*b_strip**4*d_1*i_3 \ - - 3*b_strip**6*d_y*i_4 - 56*b_strip**4*d_xy*i_5)/420/b_strip**3 + kf_mp[0, 0] = ( + ( + 5040 * d_x * i_1 + - 504 * b_strip**2 * d_1 * i_2 + - 504 * b_strip**2 * d_1 * i_3 + + 156 * b_strip**4 * d_y * i_4 + + 2016 * b_strip**2 * d_xy * i_5 + ) + / 420 + / b_strip**3 + ) + kf_mp[0, 1] = ( + ( + 2520 * b_strip * d_x * i_1 + - 462 * b_strip**3 * d_1 * i_2 + - 42 * b_strip**3 * d_1 * i_3 + + 22 * b_strip**5 * d_y * i_4 + + 168 * b_strip**3 * d_xy * i_5 + ) + / 420 + / b_strip**3 + ) + kf_mp[0, 2] = ( + ( + -5040 * d_x * i_1 + + 504 * b_strip**2 * d_1 * i_2 + + 504 * b_strip**2 * d_1 * i_3 + + 54 * b_strip**4 * d_y * i_4 + - 2016 * b_strip**2 * d_xy * i_5 + ) + / 420 + / b_strip**3 + ) + kf_mp[0, 3] = ( + ( + 2520 * b_strip * d_x * i_1 + - 42 * b_strip**3 * d_1 * i_2 + - 42 * b_strip**3 * d_1 * i_3 + - 13 * b_strip**5 * d_y * i_4 + + 168 * b_strip**3 * d_xy * i_5 + ) + / 420 + / b_strip**3 + ) + + kf_mp[1, 0] = ( + ( + 2520 * b_strip * d_x * i_1 + - 462 * b_strip**3 * d_1 * i_3 + - 42 * b_strip**3 * d_1 * i_2 + + 22 * b_strip**5 * d_y * i_4 + + 168 * b_strip**3 * d_xy * i_5 + ) + / 420 + / b_strip**3 + ) + kf_mp[1, 1] = ( + ( + 1680 * b_strip**2 * d_x * i_1 + - 56 * b_strip**4 * d_1 * i_2 + - 56 * b_strip**4 * d_1 * i_3 + + 4 * b_strip**6 * d_y * i_4 + + 224 * b_strip**4 * d_xy * i_5 + ) + / 420 + / b_strip**3 + ) + kf_mp[1, 2] = ( + ( + -2520 * b_strip * d_x * i_1 + + 42 * b_strip**3 * d_1 * i_2 + + 42 * b_strip**3 * d_1 * i_3 + + 13 * b_strip**5 * d_y * i_4 + - 168 * b_strip**3 * d_xy * i_5 + ) + / 420 + / b_strip**3 + ) + kf_mp[1, 3] = ( + ( + 840 * b_strip**2 * d_x * i_1 + + 14 * b_strip**4 * d_1 * i_2 + + 14 * b_strip**4 * d_1 * i_3 + - 3 * b_strip**6 * d_y * i_4 + - 56 * b_strip**4 * d_xy * i_5 + ) + / 420 + / b_strip**3 + ) kf_mp[2, 0] = kf_mp[0, 2] kf_mp[2, 1] = kf_mp[1, 2] - kf_mp[2, 2] = (5040*d_x*i_1 - 504*b_strip**2*d_1*i_2 - 504*b_strip**2*d_1*i_3 \ - + 156*b_strip**4*d_y*i_4 + 2016*b_strip**2*d_xy*i_5)/420/b_strip**3 - kf_mp[2, 3] = (-2520*b_strip*d_x*i_1 + 462*b_strip**3*d_1*i_2 + 42*b_strip**3*d_1*i_3 \ - - 22*b_strip**5*d_y*i_4 - 168*b_strip**3*d_xy*i_5)/420/b_strip**3 + kf_mp[2, 2] = ( + ( + 5040 * d_x * i_1 + - 504 * b_strip**2 * d_1 * i_2 + - 504 * b_strip**2 * d_1 * i_3 + + 156 * b_strip**4 * d_y * i_4 + + 2016 * b_strip**2 * d_xy * i_5 + ) + / 420 + / b_strip**3 + ) + kf_mp[2, 3] = ( + ( + -2520 * b_strip * d_x * i_1 + + 462 * b_strip**3 * d_1 * i_2 + + 42 * b_strip**3 * d_1 * i_3 + - 22 * b_strip**5 * d_y * i_4 + - 168 * b_strip**3 * d_xy * i_5 + ) + / 420 + / b_strip**3 + ) kf_mp[3, 0] = kf_mp[0, 3] kf_mp[3, 1] = kf_mp[1, 3] - kf_mp[3, 2] = (-2520*b_strip*d_x*i_1 + 462*b_strip**3*d_1*i_3 + 42*b_strip**3*d_1*i_2 \ - - 22*b_strip**5*d_y*i_4 - 168*b_strip**3*d_xy*i_5)/420/b_strip**3 # not symmetric - kf_mp[3, 3] = (1680*b_strip**2*d_x*i_1 - 56*b_strip**4*d_1*i_2 - 56*b_strip**4*d_1*i_3 \ - + 4*b_strip**6*d_y*i_4 + 224*b_strip**4*d_xy*i_5)/420/b_strip**3 + kf_mp[3, 2] = ( + ( + -2520 * b_strip * d_x * i_1 + + 462 * b_strip**3 * d_1 * i_3 + + 42 * b_strip**3 * d_1 * i_2 + - 22 * b_strip**5 * d_y * i_4 + - 168 * b_strip**3 * d_xy * i_5 + ) + / 420 + / b_strip**3 + ) # not symmetric + kf_mp[3, 3] = ( + ( + 1680 * b_strip**2 * d_x * i_1 + - 56 * b_strip**4 * d_1 * i_2 + - 56 * b_strip**4 * d_1 * i_3 + + 4 * b_strip**6 * d_y * i_4 + + 224 * b_strip**4 * d_xy * i_5 + ) + / 420 + / b_strip**3 + ) return kf_mp cdef np.ndarray calc_gm_mp( - double u_i, double u_j, double b_strip, double length, double ty_1, double ty_2, - double i_4, double i_5 + double u_i, + double u_j, + double b_strip, + double length, + double ty_1, + double ty_2, + double i_4, + double i_5, ): - """Calculate the membrane geometric stiffness sub-matrix, used in the assembly of local + """Calculate the membrane geometric stiffness sub-matrix, used in the assembly of local geometric stiffness matrices Args: @@ -642,14 +803,14 @@ cdef np.ndarray calc_gm_mp( cdef np.ndarray[np.double_t, ndim=2] gm_mp = np.zeros((4, 4)) # assemble the matrix of gm_mp (symmetric membrane stability matrix) - gm_mp[0, 0] = b_strip * (3*ty_1 + ty_2) * i_5 / 12 - gm_mp[0, 2] = b_strip * (ty_1+ty_2) * i_5 / 12 + gm_mp[0, 0] = b_strip * (3 * ty_1 + ty_2) * i_5 / 12 + gm_mp[0, 2] = b_strip * (ty_1 + ty_2) * i_5 / 12 gm_mp[2, 0] = gm_mp[0, 2] - gm_mp[1, 1] = b_strip * length**2 * (3*ty_1 + ty_2) * i_4 / 12 / u_i / u_j - gm_mp[1, 3] = b_strip * length**2 * (ty_1+ty_2) * i_4 / 12 / u_i / u_j + gm_mp[1, 1] = b_strip * length**2 * (3 * ty_1 + ty_2) * i_4 / 12 / u_i / u_j + gm_mp[1, 3] = b_strip * length**2 * (ty_1 + ty_2) * i_4 / 12 / u_i / u_j gm_mp[3, 1] = gm_mp[1, 3] - gm_mp[2, 2] = b_strip * (ty_1 + 3*ty_2) * i_5 / 12 - gm_mp[3, 3] = b_strip * length**2 * (ty_1 + 3*ty_2) * i_4 / 12 / u_i / u_j + gm_mp[2, 2] = b_strip * (ty_1 + 3 * ty_2) * i_5 / 12 + gm_mp[3, 3] = b_strip * length**2 * (ty_1 + 3 * ty_2) * i_4 / 12 / u_i / u_j return gm_mp @@ -675,27 +836,27 @@ cdef np.ndarray calc_gf_mp(double ty_1, double ty_2, double b_strip, double i_5) cdef np.ndarray[np.double_t, ndim=2] gf_mp = np.zeros((4, 4)) # assemble the matrix of gf_mp (symmetric flexural stability matrix) - gf_mp[0, 0] = (10*ty_1 + 3*ty_2) * b_strip * i_5 / 35 - gf_mp[0, 1] = (15*ty_1 + 7*ty_2) * b_strip**2 * i_5 / 210 / 2 + gf_mp[0, 0] = (10 * ty_1 + 3 * ty_2) * b_strip * i_5 / 35 + gf_mp[0, 1] = (15 * ty_1 + 7 * ty_2) * b_strip**2 * i_5 / 210 / 2 gf_mp[1, 0] = gf_mp[0, 1] - gf_mp[0, 2] = 9 * (ty_1+ty_2) * b_strip * i_5 / 140 + gf_mp[0, 2] = 9 * (ty_1 + ty_2) * b_strip * i_5 / 140 gf_mp[2, 0] = gf_mp[0, 2] - gf_mp[0, 3] = -(7*ty_1 + 6*ty_2) * b_strip**2 * i_5 / 420 + gf_mp[0, 3] = -(7 * ty_1 + 6 * ty_2) * b_strip**2 * i_5 / 420 gf_mp[3, 0] = gf_mp[0, 3] - gf_mp[1, 1] = (5*ty_1 + 3*ty_2) * b_strip**3 * i_5 / 2 / 420 - gf_mp[1, 2] = (6*ty_1 + 7*ty_2) * b_strip**2 * i_5 / 420 + gf_mp[1, 1] = (5 * ty_1 + 3 * ty_2) * b_strip**3 * i_5 / 2 / 420 + gf_mp[1, 2] = (6 * ty_1 + 7 * ty_2) * b_strip**2 * i_5 / 420 gf_mp[2, 1] = gf_mp[1, 2] gf_mp[1, 3] = -(ty_1 + ty_2) * b_strip**3 * i_5 / 140 / 2 gf_mp[3, 1] = gf_mp[1, 3] - gf_mp[2, 2] = (3*ty_1 + 10*ty_2) * b_strip * i_5 / 35 - gf_mp[2, 3] = -(7*ty_1 + 15*ty_2) * b_strip**2 * i_5 / 420 + gf_mp[2, 2] = (3 * ty_1 + 10 * ty_2) * b_strip * i_5 / 35 + gf_mp[2, 3] = -(7 * ty_1 + 15 * ty_2) * b_strip**2 * i_5 / 420 gf_mp[3, 2] = gf_mp[2, 3] - gf_mp[3, 3] = (3*ty_1 + 5*ty_2) * b_strip**3 * i_5 / 420 / 2 + gf_mp[3, 3] = (3 * ty_1 + 5 * ty_2) * b_strip**3 * i_5 / 420 / 2 return gf_mp -cdef bc_i1_5(str b_c, double m_i, double m_j, double length): +cdef tuple bc_i1_5(str b_c, double m_i, double m_j, double length): """Calculate the 5 undetermined parameters i_1,i_2,i_3,i_4,i_5 for local elastic and geometric stiffness matrices. @@ -719,22 +880,22 @@ cdef bc_i1_5(str b_c, double m_i, double m_j, double length): calculation of i_4 is the integration of y_m''*Yn'' from 0 to length calculation of i_5 is the integration of y_m'*Yn' from 0 to length """ - cdef double i_1 = 0 - cdef double i_2 = 0 - cdef double i_3 = 0 - cdef double i_4 = 0 - cdef double i_5 = 0 + cdef double i_1 = double(0.0) + cdef double i_2 = double(0.0) + cdef double i_3 = double(0.0) + cdef double i_4 = double(0.0) + cdef double i_5 = double(0.0) - if b_c == 'S-S': + if b_c == "S-S": # For simply-pimply supported boundary condition at loaded edges if m_i == m_j: i_1 = length / 2 - i_2 = -m_i**2 * np.pi**2 / length / 2 - i_3 = -m_j**2 * np.pi**2 / length / 2 + i_2 = -(m_i**2) * np.pi**2 / length / 2 + i_3 = -(m_j**2) * np.pi**2 / length / 2 i_4 = np.pi**4 * m_i**4 / 2 / length**3 i_5 = np.pi**2 * m_i**2 / 2 / length - elif b_c == 'C-C': + elif b_c == "C-C": # For Clamped-clamped boundary condition at loaded edges # calculation of i_1 is the integration of y_m*Yn from 0 to length if m_i == m_j: @@ -744,61 +905,64 @@ cdef bc_i1_5(str b_c, double m_i, double m_j, double length): i_1 = length / 4 i_2 = -(m_i**2 + 1) * np.pi**2 / 4 / length i_3 = -(m_j**2 + 1) * np.pi**2 / 4 / length - i_4 = np.pi**4 * ((m_i**2 + 1)**2 + 4 * m_i**2) / 4 / length**3 + i_4 = np.pi**4 * ((m_i**2 + 1) ** 2 + 4 * m_i**2) / 4 / length**3 i_5 = (1 + m_i**2) * np.pi**2 / 4 / length else: if m_i - m_j == 2: i_1 = -length / 8 i_2 = (m_i**2 + 1) * np.pi**2 / 8 / length - m_i * np.pi**2 / 4 / length i_3 = (m_j**2 + 1) * np.pi**2 / 8 / length + m_j * np.pi**2 / 4 / length - i_4 = -(m_i - 1)**2 * (m_j + 1)**2 * np.pi**4 / 8 / length**3 - i_5 = -(1 + m_i*m_j) * np.pi**2 / 8 / length + i_4 = -((m_i - 1) ** 2) * (m_j + 1) ** 2 * np.pi**4 / 8 / length**3 + i_5 = -(1 + m_i * m_j) * np.pi**2 / 8 / length elif m_i - m_j == -2: i_1 = -length / 8 i_2 = (m_i**2 + 1) * np.pi**2 / 8 / length + m_i * np.pi**2 / 4 / length i_3 = (m_j**2 + 1) * np.pi**2 / 8 / length - m_j * np.pi**2 / 4 / length - i_4 = -(m_i + 1)**2 * (m_j - 1)**2 * np.pi**4 / 8 / length**3 - i_5 = -(1 + m_i*m_j) * np.pi**2 / 8 / length + i_4 = -((m_i + 1) ** 2) * (m_j - 1) ** 2 * np.pi**4 / 8 / length**3 + i_5 = -(1 + m_i * m_j) * np.pi**2 / 8 / length - elif b_c == 'S-C' or b_c == 'C-S': + elif b_c in ("S-C", "C-S"): # For simply-clamped supported boundary condition at loaded edges # calculation of i_1 is the integration of y_m*Yn from 0 to length if m_i == m_j: - i_1 = (1 + (m_i + 1)**2 / m_i**2) * length / 2 - i_2 = -(m_i + 1)**2 * np.pi**2 / length - i_3 = -(m_i + 1)**2 * np.pi**2 / length - i_4 = (m_i + 1)**2 * np.pi**4 * ((m_i + 1)**2 + m_i**2) / 2 / length**3 - i_5 = (1 + m_i)**2 * np.pi**2 / length + i_1 = (1 + (m_i + 1) ** 2 / m_i**2) * length / 2 + i_2 = -((m_i + 1) ** 2) * np.pi**2 / length + i_3 = -((m_i + 1) ** 2) * np.pi**2 / length + i_4 = (m_i + 1) ** 2 * np.pi**4 * ((m_i + 1) ** 2 + m_i**2) / 2 / length**3 + i_5 = (1 + m_i) ** 2 * np.pi**2 / length else: if m_i - m_j == 1: - i_1 = (m_i+1) * length / 2 / m_i + i_1 = (m_i + 1) * length / 2 / m_i i_2 = -(m_i + 1) * m_i * np.pi**2 / 2 / length - i_3 = -(m_j + 1)**2 * np.pi**2 * (m_i+1) / 2 / length / m_i - i_4 = (m_i+1) * m_i * (m_j + 1)**2 * np.pi**4 / 2 / length**3 - i_5 = (1+m_i) * (1+m_j) * np.pi**2 / 2 / length + i_3 = -((m_j + 1) ** 2) * np.pi**2 * (m_i + 1) / 2 / length / m_i + i_4 = (m_i + 1) * m_i * (m_j + 1) ** 2 * np.pi**4 / 2 / length**3 + i_5 = (1 + m_i) * (1 + m_j) * np.pi**2 / 2 / length elif m_i - m_j == -1: - i_1 = (m_j+1) * length / 2 / m_j - i_2 = -(m_i + 1)**2 * np.pi**2 * (m_j+1) / 2 / length / m_j + i_1 = (m_j + 1) * length / 2 / m_j + i_2 = -((m_i + 1) ** 2) * np.pi**2 * (m_j + 1) / 2 / length / m_j i_3 = -(m_j + 1) * m_j * np.pi**2 / 2 / length - i_4 = (m_i + 1)**2 * m_j * (m_j+1) * np.pi**4 / 2 / length**3 - i_5 = (1+m_i) * (1+m_j) * np.pi**2 / 2 / length + i_4 = (m_i + 1) ** 2 * m_j * (m_j + 1) * np.pi**4 / 2 / length**3 + i_5 = (1 + m_i) * (1 + m_j) * np.pi**2 / 2 / length - elif b_c == 'C-F' or b_c == 'F-C': + elif b_c in ("C-F", "F-C"): # For clamped-free supported boundary condition at loaded edges # calculation of i_1 is the integration of y_m*Yn from 0 to length if m_i == m_j: - i_1 = 3*length/2 - 2 * length * (-1)**(m_i - 1) / (m_i - 1/2) / np.pi - i_2 = (m_i - 1/2)**2 * np.pi**2 * ((-1)**(m_i - 1) / (m_i - 1/2) / np.pi - 1/2) / length - i_3 = (m_j - 1/2)**2 * np.pi**2 * ((-1)**(m_j - 1) / (m_j - 1/2) / np.pi - 1/2) / length - i_4 = (m_i - 1/2)**4 * np.pi**4 / 2 / length**3 - i_5 = (m_i - 1/2)**2 * np.pi**2 / 2 / length + i_1 = 3 * length / 2 - 2 * length * (-1) ** (m_i - 1) / (m_i - 1 / 2) / np.pi + i_2 = (m_i - 1 / 2) ** 2 * np.pi**2 * ((-1) ** (m_i - 1) / (m_i - 1 / 2) / np.pi - 1 / 2) / length + i_3 = (m_j - 1 / 2) ** 2 * np.pi**2 * ((-1) ** (m_j - 1) / (m_j - 1 / 2) / np.pi - 1 / 2) / length + i_4 = (m_i - 1 / 2) ** 4 * np.pi**4 / 2 / length**3 + i_5 = (m_i - 1 / 2) ** 2 * np.pi**2 / 2 / length else: - i_1 = length - length*(-1)**(m_i - 1) / (m_i - 1/2) / np.pi \ - - length*(-1)**(m_j - 1) / (m_j - 1/2) / np.pi - i_2 = (m_i - 1/2)**2 * np.pi**2 * ((-1)**(m_i - 1) / (m_i - 1/2) / np.pi) / length - i_3 = (m_j - 1/2)**2 * np.pi**2 * ((-1)**(m_j - 1) / (m_j - 1/2) / np.pi) / length + i_1 = ( + length + - length * (-1) ** (m_i - 1) / (m_i - 1 / 2) / np.pi + - length * (-1) ** (m_j - 1) / (m_j - 1 / 2) / np.pi + ) + i_2 = (m_i - 1 / 2) ** 2 * np.pi**2 * ((-1) ** (m_i - 1) / (m_i - 1 / 2) / np.pi) / length + i_3 = (m_j - 1 / 2) ** 2 * np.pi**2 * ((-1) ** (m_j - 1) / (m_j - 1 / 2) / np.pi) / length - elif b_c == 'C-G' or b_c == 'G-C': + elif b_c in ("C-G", "G-C"): # For clamped-guided supported boundary condition at loaded edges # calculation of i_1 is the integration of y_m*Yn from 0 to length if m_i == m_j: @@ -806,30 +970,27 @@ cdef bc_i1_5(str b_c, double m_i, double m_j, double length): i_1 = 3 * length / 8 else: i_1 = length / 4 - i_2 = -((m_i - 1/2)**2 + 1/4) * np.pi**2 / length / 4 - i_3 = -((m_i - 1/2)**2 + 1/4) * np.pi**2 / length / 4 - i_4 = ((m_i - 1/2)**2 - + 1/4)**2 * np.pi**4 / 4 / length**3 + (m_i - 1/2)**2 * np.pi**4 / 4 / length**3 - i_5 = (m_i - 1/2)**2 * np.pi**2 / length / 4 + np.pi**2 / 16 / length + i_2 = -((m_i - 1 / 2) ** 2 + 1 / 4) * np.pi**2 / length / 4 + i_3 = -((m_i - 1 / 2) ** 2 + 1 / 4) * np.pi**2 / length / 4 + i_4 = ((m_i - 1 / 2) ** 2 + 1 / 4) ** 2 * np.pi**4 / 4 / length**3 + ( + m_i - 1 / 2 + ) ** 2 * np.pi**4 / 4 / length**3 + i_5 = (m_i - 1 / 2) ** 2 * np.pi**2 / length / 4 + np.pi**2 / 16 / length else: if m_i - m_j == 1: i_1 = -length / 8 - i_2 = ((m_i - 1/2)**2 - + 1/4) * np.pi**2 / length / 8 - (m_i - 1/2) * np.pi**2 / length / 8 - i_3 = ((m_j - 1/2)**2 - + 1/4) * np.pi**2 / length / 8 + (m_j - 1/2) * np.pi**2 / length / 8 - i_4 = -m_j**4 * np.pi**4 / 8 / length**3 - i_5 = -m_j**2 * np.pi**2 / 8 / length + i_2 = ((m_i - 1 / 2) ** 2 + 1 / 4) * np.pi**2 / length / 8 - (m_i - 1 / 2) * np.pi**2 / length / 8 + i_3 = ((m_j - 1 / 2) ** 2 + 1 / 4) * np.pi**2 / length / 8 + (m_j - 1 / 2) * np.pi**2 / length / 8 + i_4 = -(m_j**4) * np.pi**4 / 8 / length**3 + i_5 = -(m_j**2) * np.pi**2 / 8 / length elif m_i - m_j == -1: i_1 = -length / 8 - i_2 = ((m_i - 1/2)**2 - + 1/4) * np.pi**2 / length / 8 + (m_i - 1/2) * np.pi**2 / length / 8 - i_3 = ((m_j - 1/2)**2 - + 1/4) * np.pi**2 / length / 8 - (m_j - 1/2) * np.pi**2 / length / 8 - i_4 = -m_i**4 * np.pi**4 / 8 / length**3 - i_5 = -m_i**2 * np.pi**2 / 8 / length + i_2 = ((m_i - 1 / 2) ** 2 + 1 / 4) * np.pi**2 / length / 8 + (m_i - 1 / 2) * np.pi**2 / length / 8 + i_3 = ((m_j - 1 / 2) ** 2 + 1 / 4) * np.pi**2 / length / 8 - (m_j - 1 / 2) * np.pi**2 / length / 8 + i_4 = -(m_i**4) * np.pi**4 / 8 / length**3 + i_5 = -(m_i**2) * np.pi**2 / 8 / length - return [i_1, i_2, i_3, i_4, i_5] + return (i_1, i_2, i_3, i_4, i_5) cpdef np.ndarray trans(float alpha, int total_m): @@ -841,14 +1002,22 @@ cpdef np.ndarray trans(float alpha, int total_m): Returns: gamma (np.ndarray): transformation matrix - + Zhanjie 2008 modified by Z. Li, Aug. 09, 2009 """ - cdef np.ndarray[np.double_t, ndim=2] gam = np.array([[np.cos(alpha), 0, 0, 0, -np.sin(alpha), 0, 0, 0], [0, 1, 0, 0, 0, 0, 0, 0], - [0, 0, np.cos(alpha), 0, 0, 0, -np.sin(alpha), 0], [0, 0, 0, 1, 0, 0, 0, 0], - [np.sin(alpha), 0, 0, 0, np.cos(alpha), 0, 0, 0], [0, 0, 0, 0, 0, 1, 0, 0], - [0, 0, np.sin(alpha), 0, 0, 0, np.cos(alpha), 0], [0, 0, 0, 0, 0, 0, 0, 1]], dtype=np.double) + cdef np.ndarray[np.double_t, ndim=2] gam = np.array( + [ + [np.cos(alpha), 0, 0, 0, -np.sin(alpha), 0, 0, 0], + [0, 1, 0, 0, 0, 0, 0, 0], + [0, 0, np.cos(alpha), 0, 0, 0, -np.sin(alpha), 0], + [0, 0, 0, 1, 0, 0, 0, 0], + [np.sin(alpha), 0, 0, 0, np.cos(alpha), 0, 0, 0], + [0, 0, 0, 0, 0, 1, 0, 0], + [0, 0, np.sin(alpha), 0, 0, 0, np.cos(alpha), 0], + [0, 0, 0, 0, 0, 0, 0, 1], + ] + ) if total_m == 1: return gam @@ -857,14 +1026,20 @@ cpdef np.ndarray trans(float alpha, int total_m): # extend to multi-m cdef int i for i in range(0, total_m): - gamma[8 * i:8 * (i+1), 8 * i:8 * (i+1)] = gam + gamma[8 * i : 8 * (i + 1), 8 * i : 8 * (i + 1)] = gam return gamma -cdef assemble( - np.ndarray k_global, np.ndarray kg_global, np.ndarray k_local, np.ndarray kg_local, - int node_i, int node_j, int n_nodes, np.ndarray m_a +cdef tuple assemble( + np.ndarray k_global, + np.ndarray kg_global, + np.ndarray k_local, + np.ndarray kg_local, + int node_i, + int node_j, + int n_nodes, + np.ndarray m_a, ): """Add the element contribution to the global stiffness matrix @@ -875,7 +1050,7 @@ cdef assemble( one used in original CUFSM for single longitudinal term m in the DOF order [u1 v1...un vn w1 01...wn 0n]m'. kg_local (np.ndarray): local geometric stiffness matrix. Each submatrix is similar to the - one used in original CUFSM for single longitudinal term m in the DOF order + one used in original CUFSM for single longitudinal term m in the DOF order [u1 v1...un vn w1 01...wn 0n]m'. node_i (int): node number node_j (int): node number @@ -932,122 +1107,186 @@ cdef assemble( for i in range(0, total_m): for j in range(0, total_m): # Submatrices for the initial stiffness - k11 = k_local[8 * i:8*i + 2, 8 * j:8*j + 2] - k12 = k_local[8 * i:8*i + 2, 8*j + 2:8*j + 4] - k13 = k_local[8 * i:8*i + 2, 8*j + 4:8*j + 6] - k14 = k_local[8 * i:8*i + 2, 8*j + 6:8*j + 8] - k21 = k_local[8*i + 2:8*i + 4, 8 * j:8*j + 2] - k22 = k_local[8*i + 2:8*i + 4, 8*j + 2:8*j + 4] - k23 = k_local[8*i + 2:8*i + 4, 8*j + 4:8*j + 6] - k24 = k_local[8*i + 2:8*i + 4, 8*j + 6:8*j + 8] - k31 = k_local[8*i + 4:8*i + 6, 8 * j:8*j + 2] - k32 = k_local[8*i + 4:8*i + 6, 8*j + 2:8*j + 4] - k33 = k_local[8*i + 4:8*i + 6, 8*j + 4:8*j + 6] - k34 = k_local[8*i + 4:8*i + 6, 8*j + 6:8*j + 8] - k41 = k_local[8*i + 6:8*i + 8, 8 * j:8*j + 2] - k42 = k_local[8*i + 6:8*i + 8, 8*j + 2:8*j + 4] - k43 = k_local[8*i + 6:8*i + 8, 8*j + 4:8*j + 6] - k44 = k_local[8*i + 6:8*i + 8, 8*j + 6:8*j + 8] - - k_global[4*n_nodes*i+(node_i+1)*2-2:4*n_nodes*i+(node_i+1)*2, \ - 4*n_nodes*j+(node_i+1)*2-2:4*n_nodes*j+(node_i+1)*2] += k11 - k_global[4*n_nodes*i+(node_i+1)*2-2:4*n_nodes*i+(node_i+1)*2, \ - 4*n_nodes*j+(node_j+1)*2-2:4*n_nodes*j+(node_j+1)*2] += k12 - k_global[4*n_nodes*i+(node_j+1)*2-2:4*n_nodes*i+(node_j+1)*2, \ - 4*n_nodes*j+(node_i+1)*2-2:4*n_nodes*j+(node_i+1)*2] += k21 - k_global[4*n_nodes*i+(node_j+1)*2-2:4*n_nodes*i+(node_j+1)*2, \ - 4*n_nodes*j+(node_j+1)*2-2:4*n_nodes*j+(node_j+1)*2] += k22 - - k_global[4*n_nodes*i+skip+(node_i+1)*2-2:4*n_nodes*i+skip+(node_i+1)*2, \ - 4*n_nodes*j+skip+(node_i+1)*2-2:4*n_nodes*j+skip+(node_i+1)*2] += k33 - k_global[4*n_nodes*i+skip+(node_i+1)*2-2:4*n_nodes*i+skip+(node_i+1)*2, \ - 4*n_nodes*j+skip+(node_j+1)*2-2:4*n_nodes*j+skip+(node_j+1)*2] += k34 - k_global[4*n_nodes*i+skip+(node_j+1)*2-2:4*n_nodes*i+skip+(node_j+1)*2, \ - 4*n_nodes*j+skip+(node_i+1)*2-2:4*n_nodes*j+skip+(node_i+1)*2] += k43 - k_global[4*n_nodes*i+skip+(node_j+1)*2-2:4*n_nodes*i+skip+(node_j+1)*2, \ - 4*n_nodes*j+skip+(node_j+1)*2-2:4*n_nodes*j+skip+(node_j+1)*2] += k44 - - k_global[4*n_nodes*i+(node_i+1)*2-2:4*n_nodes*i+(node_i+1)*2, \ - 4*n_nodes*j+skip+(node_i+1)*2-2:4*n_nodes*j+skip+(node_i+1)*2] += k13 - k_global[4*n_nodes*i+(node_i+1)*2-2:4*n_nodes*i+(node_i+1)*2, \ - 4*n_nodes*j+skip+(node_j+1)*2-2:4*n_nodes*j+skip+(node_j+1)*2] += k14 - k_global[4*n_nodes*i+(node_j+1)*2-2:4*n_nodes*i+(node_j+1)*2, \ - 4*n_nodes*j+skip+(node_i+1)*2-2:4*n_nodes*j+skip+(node_i+1)*2] += k23 - k_global[4*n_nodes*i+(node_j+1)*2-2:4*n_nodes*i+(node_j+1)*2, \ - 4*n_nodes*j+skip+(node_j+1)*2-2:4*n_nodes*j+skip+(node_j+1)*2] += k24 - - k_global[4*n_nodes*i+skip+(node_i+1)*2-2:4*n_nodes*i+skip+(node_i+1)*2, \ - 4*n_nodes*j+(node_i+1)*2-2:4*n_nodes*j+(node_i+1)*2] += k31 - k_global[4*n_nodes*i+skip+(node_i+1)*2-2:4*n_nodes*i+skip+(node_i+1)*2, \ - 4*n_nodes*j+(node_j+1)*2-2:4*n_nodes*j+(node_j+1)*2] += k32 - k_global[4*n_nodes*i+skip+(node_j+1)*2-2:4*n_nodes*i+skip+(node_j+1)*2, \ - 4*n_nodes*j+(node_i+1)*2-2:4*n_nodes*j+(node_i+1)*2] += k41 - k_global[4*n_nodes*i+skip+(node_j+1)*2-2:4*n_nodes*i+skip+(node_j+1)*2, \ - 4*n_nodes*j+(node_j+1)*2-2:4*n_nodes*j+(node_j+1)*2] += k42 + k11 = k_local[8 * i : 8 * i + 2, 8 * j : 8 * j + 2] + k12 = k_local[8 * i : 8 * i + 2, 8 * j + 2 : 8 * j + 4] + k13 = k_local[8 * i : 8 * i + 2, 8 * j + 4 : 8 * j + 6] + k14 = k_local[8 * i : 8 * i + 2, 8 * j + 6 : 8 * j + 8] + k21 = k_local[8 * i + 2 : 8 * i + 4, 8 * j : 8 * j + 2] + k22 = k_local[8 * i + 2 : 8 * i + 4, 8 * j + 2 : 8 * j + 4] + k23 = k_local[8 * i + 2 : 8 * i + 4, 8 * j + 4 : 8 * j + 6] + k24 = k_local[8 * i + 2 : 8 * i + 4, 8 * j + 6 : 8 * j + 8] + k31 = k_local[8 * i + 4 : 8 * i + 6, 8 * j : 8 * j + 2] + k32 = k_local[8 * i + 4 : 8 * i + 6, 8 * j + 2 : 8 * j + 4] + k33 = k_local[8 * i + 4 : 8 * i + 6, 8 * j + 4 : 8 * j + 6] + k34 = k_local[8 * i + 4 : 8 * i + 6, 8 * j + 6 : 8 * j + 8] + k41 = k_local[8 * i + 6 : 8 * i + 8, 8 * j : 8 * j + 2] + k42 = k_local[8 * i + 6 : 8 * i + 8, 8 * j + 2 : 8 * j + 4] + k43 = k_local[8 * i + 6 : 8 * i + 8, 8 * j + 4 : 8 * j + 6] + k44 = k_local[8 * i + 6 : 8 * i + 8, 8 * j + 6 : 8 * j + 8] + + k_global[ + 4 * n_nodes * i + (node_i + 1) * 2 - 2 : 4 * n_nodes * i + (node_i + 1) * 2, + 4 * n_nodes * j + (node_i + 1) * 2 - 2 : 4 * n_nodes * j + (node_i + 1) * 2, + ] += k11 + k_global[ + 4 * n_nodes * i + (node_i + 1) * 2 - 2 : 4 * n_nodes * i + (node_i + 1) * 2, + 4 * n_nodes * j + (node_j + 1) * 2 - 2 : 4 * n_nodes * j + (node_j + 1) * 2, + ] += k12 + k_global[ + 4 * n_nodes * i + (node_j + 1) * 2 - 2 : 4 * n_nodes * i + (node_j + 1) * 2, + 4 * n_nodes * j + (node_i + 1) * 2 - 2 : 4 * n_nodes * j + (node_i + 1) * 2, + ] += k21 + k_global[ + 4 * n_nodes * i + (node_j + 1) * 2 - 2 : 4 * n_nodes * i + (node_j + 1) * 2, + 4 * n_nodes * j + (node_j + 1) * 2 - 2 : 4 * n_nodes * j + (node_j + 1) * 2, + ] += k22 + + k_global[ + 4 * n_nodes * i + skip + (node_i + 1) * 2 - 2 : 4 * n_nodes * i + skip + (node_i + 1) * 2, + 4 * n_nodes * j + skip + (node_i + 1) * 2 - 2 : 4 * n_nodes * j + skip + (node_i + 1) * 2, + ] += k33 + k_global[ + 4 * n_nodes * i + skip + (node_i + 1) * 2 - 2 : 4 * n_nodes * i + skip + (node_i + 1) * 2, + 4 * n_nodes * j + skip + (node_j + 1) * 2 - 2 : 4 * n_nodes * j + skip + (node_j + 1) * 2, + ] += k34 + k_global[ + 4 * n_nodes * i + skip + (node_j + 1) * 2 - 2 : 4 * n_nodes * i + skip + (node_j + 1) * 2, + 4 * n_nodes * j + skip + (node_i + 1) * 2 - 2 : 4 * n_nodes * j + skip + (node_i + 1) * 2, + ] += k43 + k_global[ + 4 * n_nodes * i + skip + (node_j + 1) * 2 - 2 : 4 * n_nodes * i + skip + (node_j + 1) * 2, + 4 * n_nodes * j + skip + (node_j + 1) * 2 - 2 : 4 * n_nodes * j + skip + (node_j + 1) * 2, + ] += k44 + + k_global[ + 4 * n_nodes * i + (node_i + 1) * 2 - 2 : 4 * n_nodes * i + (node_i + 1) * 2, + 4 * n_nodes * j + skip + (node_i + 1) * 2 - 2 : 4 * n_nodes * j + skip + (node_i + 1) * 2, + ] += k13 + k_global[ + 4 * n_nodes * i + (node_i + 1) * 2 - 2 : 4 * n_nodes * i + (node_i + 1) * 2, + 4 * n_nodes * j + skip + (node_j + 1) * 2 - 2 : 4 * n_nodes * j + skip + (node_j + 1) * 2, + ] += k14 + k_global[ + 4 * n_nodes * i + (node_j + 1) * 2 - 2 : 4 * n_nodes * i + (node_j + 1) * 2, + 4 * n_nodes * j + skip + (node_i + 1) * 2 - 2 : 4 * n_nodes * j + skip + (node_i + 1) * 2, + ] += k23 + k_global[ + 4 * n_nodes * i + (node_j + 1) * 2 - 2 : 4 * n_nodes * i + (node_j + 1) * 2, + 4 * n_nodes * j + skip + (node_j + 1) * 2 - 2 : 4 * n_nodes * j + skip + (node_j + 1) * 2, + ] += k24 + + k_global[ + 4 * n_nodes * i + skip + (node_i + 1) * 2 - 2 : 4 * n_nodes * i + skip + (node_i + 1) * 2, + 4 * n_nodes * j + (node_i + 1) * 2 - 2 : 4 * n_nodes * j + (node_i + 1) * 2, + ] += k31 + k_global[ + 4 * n_nodes * i + skip + (node_i + 1) * 2 - 2 : 4 * n_nodes * i + skip + (node_i + 1) * 2, + 4 * n_nodes * j + (node_j + 1) * 2 - 2 : 4 * n_nodes * j + (node_j + 1) * 2, + ] += k32 + k_global[ + 4 * n_nodes * i + skip + (node_j + 1) * 2 - 2 : 4 * n_nodes * i + skip + (node_j + 1) * 2, + 4 * n_nodes * j + (node_i + 1) * 2 - 2 : 4 * n_nodes * j + (node_i + 1) * 2, + ] += k41 + k_global[ + 4 * n_nodes * i + skip + (node_j + 1) * 2 - 2 : 4 * n_nodes * i + skip + (node_j + 1) * 2, + 4 * n_nodes * j + (node_j + 1) * 2 - 2 : 4 * n_nodes * j + (node_j + 1) * 2, + ] += k42 # Submatrices for the initial stiffness - kg11 = kg_local[8 * i:8*i + 2, 8 * j:8*j + 2] - kg12 = kg_local[8 * i:8*i + 2, 8*j + 2:8*j + 4] - kg13 = kg_local[8 * i:8*i + 2, 8*j + 4:8*j + 6] - kg14 = kg_local[8 * i:8*i + 2, 8*j + 6:8*j + 8] - kg21 = kg_local[8*i + 2:8*i + 4, 8 * j:8*j + 2] - kg22 = kg_local[8*i + 2:8*i + 4, 8*j + 2:8*j + 4] - kg23 = kg_local[8*i + 2:8*i + 4, 8*j + 4:8*j + 6] - kg24 = kg_local[8*i + 2:8*i + 4, 8*j + 6:8*j + 8] - kg31 = kg_local[8*i + 4:8*i + 6, 8 * j:8*j + 2] - kg32 = kg_local[8*i + 4:8*i + 6, 8*j + 2:8*j + 4] - kg33 = kg_local[8*i + 4:8*i + 6, 8*j + 4:8*j + 6] - kg34 = kg_local[8*i + 4:8*i + 6, 8*j + 6:8*j + 8] - kg41 = kg_local[8*i + 6:8*i + 8, 8 * j:8*j + 2] - kg42 = kg_local[8*i + 6:8*i + 8, 8*j + 2:8*j + 4] - kg43 = kg_local[8*i + 6:8*i + 8, 8*j + 4:8*j + 6] - kg44 = kg_local[8*i + 6:8*i + 8, 8*j + 6:8*j + 8] - - kg_global[4*n_nodes*i + (node_i+1) * 2 - 2:4*n_nodes*i + (node_i+1) * 2, - 4*n_nodes*j + (node_i+1) * 2 - 2:4*n_nodes*j + (node_i+1) * 2] += kg11 - kg_global[4*n_nodes*i + (node_i+1) * 2 - 2:4*n_nodes*i + (node_i+1) * 2, - 4*n_nodes*j + (node_j+1) * 2 - 2:4*n_nodes*j + (node_j+1) * 2] += kg12 - kg_global[4*n_nodes*i + (node_j+1) * 2 - 2:4*n_nodes*i + (node_j+1) * 2, - 4*n_nodes*j + (node_i+1) * 2 - 2:4*n_nodes*j + (node_i+1) * 2] += kg21 - kg_global[4*n_nodes*i + (node_j+1) * 2 - 2:4*n_nodes*i + (node_j+1) * 2, - 4*n_nodes*j + (node_j+1) * 2 - 2:4*n_nodes*j + (node_j+1) * 2] += kg22 - - kg_global[4*n_nodes*i + skip + (node_i+1) * 2 - 2:4*n_nodes*i + skip + (node_i+1) * 2, - 4*n_nodes*j + skip + (node_i+1) * 2 - 2:4*n_nodes*j + skip - + (node_i+1) * 2] += kg33 - kg_global[4*n_nodes*i + skip + (node_i+1) * 2 - 2:4*n_nodes*i + skip + (node_i+1) * 2, - 4*n_nodes*j + skip + (node_j+1) * 2 - 2:4*n_nodes*j + skip - + (node_j+1) * 2] += kg34 - kg_global[4*n_nodes*i + skip + (node_j+1) * 2 - 2:4*n_nodes*i + skip + (node_j+1) * 2, - 4*n_nodes*j + skip + (node_i+1) * 2 - 2:4*n_nodes*j + skip - + (node_i+1) * 2] += kg43 - kg_global[4*n_nodes*i + skip + (node_j+1) * 2 - 2:4*n_nodes*i + skip + (node_j+1) * 2, - 4*n_nodes*j + skip + (node_j+1) * 2 - 2:4*n_nodes*j + skip - + (node_j+1) * 2] += kg44 - - kg_global[4*n_nodes*i + (node_i+1) * 2 - 2:4*n_nodes*i + (node_i+1) * 2, 4*n_nodes*j - + skip + (node_i+1) * 2 - 2:4*n_nodes*j + skip + (node_i+1) * 2] += kg13 - kg_global[4*n_nodes*i + (node_i+1) * 2 - 2:4*n_nodes*i + (node_i+1) * 2, 4*n_nodes*j - + skip + (node_j+1) * 2 - 2:4*n_nodes*j + skip + (node_j+1) * 2] += kg14 - kg_global[4*n_nodes*i + (node_j+1) * 2 - 2:4*n_nodes*i + (node_j+1) * 2, 4*n_nodes*j - + skip + (node_i+1) * 2 - 2:4*n_nodes*j + skip + (node_i+1) * 2] += kg23 - kg_global[4*n_nodes*i + (node_j+1) * 2 - 2:4*n_nodes*i + (node_j+1) * 2, 4*n_nodes*j - + skip + (node_j+1) * 2 - 2:4*n_nodes*j + skip + (node_j+1) * 2] += kg24 - - kg_global[4*n_nodes*i + skip + (node_i+1) * 2 - 2:4*n_nodes*i + skip + (node_i+1) * 2, - 4*n_nodes*j + (node_i+1) * 2 - 2:4*n_nodes*j + (node_i+1) * 2] += kg31 - kg_global[4*n_nodes*i + skip + (node_i+1) * 2 - 2:4*n_nodes*i + skip + (node_i+1) * 2, - 4*n_nodes*j + (node_j+1) * 2 - 2:4*n_nodes*j + (node_j+1) * 2] += kg32 - kg_global[4*n_nodes*i + skip + (node_j+1) * 2 - 2:4*n_nodes*i + skip + (node_j+1) * 2, - 4*n_nodes*j + (node_i+1) * 2 - 2:4*n_nodes*j + (node_i+1) * 2] += kg41 - kg_global[4*n_nodes*i + skip + (node_j+1) * 2 - 2:4*n_nodes*i + skip + (node_j+1) * 2, - 4*n_nodes*j + (node_j+1) * 2 - 2:4*n_nodes*j + (node_j+1) * 2] += kg42 + kg11 = kg_local[8 * i : 8 * i + 2, 8 * j : 8 * j + 2] + kg12 = kg_local[8 * i : 8 * i + 2, 8 * j + 2 : 8 * j + 4] + kg13 = kg_local[8 * i : 8 * i + 2, 8 * j + 4 : 8 * j + 6] + kg14 = kg_local[8 * i : 8 * i + 2, 8 * j + 6 : 8 * j + 8] + kg21 = kg_local[8 * i + 2 : 8 * i + 4, 8 * j : 8 * j + 2] + kg22 = kg_local[8 * i + 2 : 8 * i + 4, 8 * j + 2 : 8 * j + 4] + kg23 = kg_local[8 * i + 2 : 8 * i + 4, 8 * j + 4 : 8 * j + 6] + kg24 = kg_local[8 * i + 2 : 8 * i + 4, 8 * j + 6 : 8 * j + 8] + kg31 = kg_local[8 * i + 4 : 8 * i + 6, 8 * j : 8 * j + 2] + kg32 = kg_local[8 * i + 4 : 8 * i + 6, 8 * j + 2 : 8 * j + 4] + kg33 = kg_local[8 * i + 4 : 8 * i + 6, 8 * j + 4 : 8 * j + 6] + kg34 = kg_local[8 * i + 4 : 8 * i + 6, 8 * j + 6 : 8 * j + 8] + kg41 = kg_local[8 * i + 6 : 8 * i + 8, 8 * j : 8 * j + 2] + kg42 = kg_local[8 * i + 6 : 8 * i + 8, 8 * j + 2 : 8 * j + 4] + kg43 = kg_local[8 * i + 6 : 8 * i + 8, 8 * j + 4 : 8 * j + 6] + kg44 = kg_local[8 * i + 6 : 8 * i + 8, 8 * j + 6 : 8 * j + 8] + + kg_global[ + 4 * n_nodes * i + (node_i + 1) * 2 - 2 : 4 * n_nodes * i + (node_i + 1) * 2, + 4 * n_nodes * j + (node_i + 1) * 2 - 2 : 4 * n_nodes * j + (node_i + 1) * 2, + ] += kg11 + kg_global[ + 4 * n_nodes * i + (node_i + 1) * 2 - 2 : 4 * n_nodes * i + (node_i + 1) * 2, + 4 * n_nodes * j + (node_j + 1) * 2 - 2 : 4 * n_nodes * j + (node_j + 1) * 2, + ] += kg12 + kg_global[ + 4 * n_nodes * i + (node_j + 1) * 2 - 2 : 4 * n_nodes * i + (node_j + 1) * 2, + 4 * n_nodes * j + (node_i + 1) * 2 - 2 : 4 * n_nodes * j + (node_i + 1) * 2, + ] += kg21 + kg_global[ + 4 * n_nodes * i + (node_j + 1) * 2 - 2 : 4 * n_nodes * i + (node_j + 1) * 2, + 4 * n_nodes * j + (node_j + 1) * 2 - 2 : 4 * n_nodes * j + (node_j + 1) * 2, + ] += kg22 + + kg_global[ + 4 * n_nodes * i + skip + (node_i + 1) * 2 - 2 : 4 * n_nodes * i + skip + (node_i + 1) * 2, + 4 * n_nodes * j + skip + (node_i + 1) * 2 - 2 : 4 * n_nodes * j + skip + (node_i + 1) * 2, + ] += kg33 + kg_global[ + 4 * n_nodes * i + skip + (node_i + 1) * 2 - 2 : 4 * n_nodes * i + skip + (node_i + 1) * 2, + 4 * n_nodes * j + skip + (node_j + 1) * 2 - 2 : 4 * n_nodes * j + skip + (node_j + 1) * 2, + ] += kg34 + kg_global[ + 4 * n_nodes * i + skip + (node_j + 1) * 2 - 2 : 4 * n_nodes * i + skip + (node_j + 1) * 2, + 4 * n_nodes * j + skip + (node_i + 1) * 2 - 2 : 4 * n_nodes * j + skip + (node_i + 1) * 2, + ] += kg43 + kg_global[ + 4 * n_nodes * i + skip + (node_j + 1) * 2 - 2 : 4 * n_nodes * i + skip + (node_j + 1) * 2, + 4 * n_nodes * j + skip + (node_j + 1) * 2 - 2 : 4 * n_nodes * j + skip + (node_j + 1) * 2, + ] += kg44 + + kg_global[ + 4 * n_nodes * i + (node_i + 1) * 2 - 2 : 4 * n_nodes * i + (node_i + 1) * 2, + 4 * n_nodes * j + skip + (node_i + 1) * 2 - 2 : 4 * n_nodes * j + skip + (node_i + 1) * 2, + ] += kg13 + kg_global[ + 4 * n_nodes * i + (node_i + 1) * 2 - 2 : 4 * n_nodes * i + (node_i + 1) * 2, + 4 * n_nodes * j + skip + (node_j + 1) * 2 - 2 : 4 * n_nodes * j + skip + (node_j + 1) * 2, + ] += kg14 + kg_global[ + 4 * n_nodes * i + (node_j + 1) * 2 - 2 : 4 * n_nodes * i + (node_j + 1) * 2, + 4 * n_nodes * j + skip + (node_i + 1) * 2 - 2 : 4 * n_nodes * j + skip + (node_i + 1) * 2, + ] += kg23 + kg_global[ + 4 * n_nodes * i + (node_j + 1) * 2 - 2 : 4 * n_nodes * i + (node_j + 1) * 2, + 4 * n_nodes * j + skip + (node_j + 1) * 2 - 2 : 4 * n_nodes * j + skip + (node_j + 1) * 2, + ] += kg24 + + kg_global[ + 4 * n_nodes * i + skip + (node_i + 1) * 2 - 2 : 4 * n_nodes * i + skip + (node_i + 1) * 2, + 4 * n_nodes * j + (node_i + 1) * 2 - 2 : 4 * n_nodes * j + (node_i + 1) * 2, + ] += kg31 + kg_global[ + 4 * n_nodes * i + skip + (node_i + 1) * 2 - 2 : 4 * n_nodes * i + skip + (node_i + 1) * 2, + 4 * n_nodes * j + (node_j + 1) * 2 - 2 : 4 * n_nodes * j + (node_j + 1) * 2, + ] += kg32 + kg_global[ + 4 * n_nodes * i + skip + (node_j + 1) * 2 - 2 : 4 * n_nodes * i + skip + (node_j + 1) * 2, + 4 * n_nodes * j + (node_i + 1) * 2 - 2 : 4 * n_nodes * j + (node_i + 1) * 2, + ] += kg41 + kg_global[ + 4 * n_nodes * i + skip + (node_j + 1) * 2 - 2 : 4 * n_nodes * i + skip + (node_j + 1) * 2, + 4 * n_nodes * j + (node_j + 1) * 2 - 2 : 4 * n_nodes * j + (node_j + 1) * 2, + ] += kg42 return k_global, kg_global cdef assemble_single( - np.ndarray k_global, np.ndarray k_local, int node_i, int node_j, int n_nodes + np.ndarray k_global, + np.ndarray k_local, + int node_i, + int node_j, + int n_nodes, ): """this routine adds the element contribution to the global stiffness matrix basically it does the same as routine 'assemble', however: @@ -1087,32 +1326,51 @@ cdef assemble_single( # the additional terms for k_global are stored in k_2_matrix cdef int skip = 2 * n_nodes - k_global[node_i * 2:node_i*2 + 2, node_i * 2:node_i*2 + 2] += k11 - k_global[node_i * 2:node_i*2 + 2, node_j * 2:node_j*2 + 2] += k12 - k_global[node_j * 2:node_j*2 + 2, node_i * 2:node_i*2 + 2] += k21 - k_global[node_j * 2:node_j*2 + 2, node_j * 2:node_j*2 + 2] += k22 - - k_global[skip + node_i*2:skip + node_i*2 + 2, skip + node_i*2:skip + node_i*2 + 2] += k33 - k_global[skip + node_i*2:skip + node_i*2 + 2, skip + node_j*2:skip + node_j*2 + 2] += k34 - k_global[skip + node_j*2:skip + node_j*2 + 2, skip + node_i*2:skip + node_i*2 + 2] += k43 - k_global[skip + node_j*2:skip + node_j*2 + 2, skip + node_j*2:skip + node_j*2 + 2] += k44 - - k_global[node_i * 2:node_i*2 + 2, skip + node_i*2:skip + node_i*2 + 2] += k13 - k_global[node_i * 2:node_i*2 + 2, skip + node_j*2:skip + node_j*2 + 2] += k14 - k_global[node_j * 2:node_j*2 + 2, skip + node_i*2:skip + node_i*2 + 2] += k23 - k_global[node_j * 2:node_j*2 + 2, skip + node_j*2:skip + node_j*2 + 2] += k24 - - k_global[skip + node_i*2:skip + node_i*2 + 2, node_i * 2:node_i*2 + 2] += k31 - k_global[skip + node_i*2:skip + node_i*2 + 2, node_j * 2:node_j*2 + 2] += k32 - k_global[skip + node_j*2:skip + node_j*2 + 2, node_i * 2:node_i*2 + 2] += k41 - k_global[skip + node_j*2:skip + node_j*2 + 2, node_j * 2:node_j*2 + 2] += k42 + k_global[node_i * 2 : node_i * 2 + 2, node_i * 2 : node_i * 2 + 2] += k11 + k_global[node_i * 2 : node_i * 2 + 2, node_j * 2 : node_j * 2 + 2] += k12 + k_global[node_j * 2 : node_j * 2 + 2, node_i * 2 : node_i * 2 + 2] += k21 + k_global[node_j * 2 : node_j * 2 + 2, node_j * 2 : node_j * 2 + 2] += k22 + + k_global[ + skip + node_i * 2 : skip + node_i * 2 + 2, + skip + node_i * 2 : skip + node_i * 2 + 2, + ] += k33 + k_global[ + skip + node_i * 2 : skip + node_i * 2 + 2, + skip + node_j * 2 : skip + node_j * 2 + 2, + ] += k34 + k_global[ + skip + node_j * 2 : skip + node_j * 2 + 2, + skip + node_i * 2 : skip + node_i * 2 + 2, + ] += k43 + k_global[ + skip + node_j * 2 : skip + node_j * 2 + 2, + skip + node_j * 2 : skip + node_j * 2 + 2, + ] += k44 + + k_global[node_i * 2 : node_i * 2 + 2, skip + node_i * 2 : skip + node_i * 2 + 2] += k13 + k_global[node_i * 2 : node_i * 2 + 2, skip + node_j * 2 : skip + node_j * 2 + 2] += k14 + k_global[node_j * 2 : node_j * 2 + 2, skip + node_i * 2 : skip + node_i * 2 + 2] += k23 + k_global[node_j * 2 : node_j * 2 + 2, skip + node_j * 2 : skip + node_j * 2 + 2] += k24 + + k_global[skip + node_i * 2 : skip + node_i * 2 + 2, node_i * 2 : node_i * 2 + 2] += k31 + k_global[skip + node_i * 2 : skip + node_i * 2 + 2, node_j * 2 : node_j * 2 + 2] += k32 + k_global[skip + node_j * 2 : skip + node_j * 2 + 2, node_i * 2 : node_i * 2 + 2] += k41 + k_global[skip + node_j * 2 : skip + node_j * 2 + 2, node_j * 2 : node_j * 2 + 2] += k42 return k_global cpdef np.ndarray spring_klocal( - double k_u, double k_v, double k_w, double k_q, double length, - str b_c, np.ndarray m_a, int discrete, double y_s + double k_u, + double k_v, + double k_w, + double k_q, + double length, + str b_c, + np.ndarray m_a, + int discrete, + double y_s, ): """Generate spring stiffness matrix (k_local) in local coordinates, modified from klocal @@ -1134,7 +1392,7 @@ cpdef np.ndarray spring_klocal( y_s (float): location of discrete spring Returns: - klocal (np.ndarray): local stiffness matrix, a total_m x total_m matrix of 8 by 8 + klocal (np.ndarray): local stiffness matrix, a total_m x total_m matrix of 8 by 8 submatrices. k_local=[k_mp]total_m x total_m block matrix each k_mp is the 8 x 8 submatrix in the DOF order [u1 v1 u2 v2 w1 theta1 w2 theta2]' @@ -1159,29 +1417,45 @@ cpdef np.ndarray spring_klocal( u_j = m_a[j] * np.pi if discrete: - [i_1, i_5] = bc_i1_5_atpoint( - b_c=b_c, m_i=m_a[i], m_j=m_a[j], length=length, y_s=y_s - ) + [i_1, i_5] = bc_i1_5_atpoint(b_c=b_c, m_i=m_a[i], m_j=m_a[j], length=length, y_s=y_s) else: # foundation spring [i_1, _, _, _, i_5] = bc_i1_5(b_c=b_c, m_i=m_a[i], m_j=m_a[j], length=length) # assemble the matrix of km_mp (membrane stiffness) km_mp = np.array( - [[k_u * i_1, 0, -k_u * i_1, 0], - [0, k_v * i_5 * length**2 / (u_i*u_j), 0, -k_v * i_5 * length**2 / (u_i*u_j)], - [-k_u * i_1, 0, k_u * i_1, 0], - [0, -k_v * i_5 * length**2 / (u_i*u_j), 0, k_v * i_5 * length**2 / (u_i*u_j)]] + [ + [k_u * i_1, 0, -k_u * i_1, 0], + [ + 0, + k_v * i_5 * length**2 / (u_i * u_j), + 0, + -k_v * i_5 * length**2 / (u_i * u_j), + ], + [-k_u * i_1, 0, k_u * i_1, 0], + [ + 0, + -k_v * i_5 * length**2 / (u_i * u_j), + 0, + k_v * i_5 * length**2 / (u_i * u_j), + ], + ] ) # assemble the matrix of kf_mp (flexural stiffness) - kf_mp = np.array([[k_w * i_1, 0, -k_w * i_1, 0], [0, k_q * i_1, 0, -k_q * i_1], - [-k_w * i_1, 0, k_w * i_1, 0], [0, -k_q * i_1, 0, k_q * i_1]]) + kf_mp = np.array( + [ + [k_w * i_1, 0, -k_w * i_1, 0], + [0, k_q * i_1, 0, -k_q * i_1], + [-k_w * i_1, 0, k_w * i_1, 0], + [0, -k_q * i_1, 0, k_q * i_1], + ] + ) - k_local[8 * i:8*i + 4, 8 * j:8*j + 4] = km_mp - k_local[8*i + 4:8 * (i+1), 8*j + 4:8 * (j+1)] = kf_mp + k_local[8 * i : 8 * i + 4, 8 * j : 8 * j + 4] = km_mp + k_local[8 * i + 4 : 8 * (i + 1), 8 * j + 4 : 8 * (j + 1)] = kf_mp return k_local -cdef bc_i1_5_atpoint(str b_c, double m_i, double m_j, double length, double y_s): +cdef tuple bc_i1_5_atpoint(str b_c, double m_i, double m_j, double length, double y_s): """Calculate the value of the longitudinal shape functions for discrete springs @@ -1211,20 +1485,24 @@ cdef bc_i1_5_atpoint(str b_c, double m_i, double m_j, double length, double y_s) cpdef np.ndarray spring_assemble( - np.ndarray k_global, np.ndarray k_local, int node_i, int node_j, int n_nodes, - np.ndarray m_a + np.ndarray k_global, + np.ndarray k_local, + int node_i, + int node_j, + int n_nodes, + np.ndarray m_a, ): """Add the (spring) contribution to the global stiffness matrix Args: - k_global (np.ndarray): global elastic stiffness matrix + k_global (np.ndarray): global elastic stiffness matrix total_m x total_m submatrices. Each submatrix is similar to the one used in original CUFSM for single longitudinal term m in the DOF order [u1 v1...un vn w1 01...wn 0n]m'. k_local (np.ndarray): local elastic stiffness matrix node_i (int): node number node_j (int): node number - n_nodes (int): total number of nodes + n_nodes (int): total number of nodes m_a (np.ndarray): number of half-wavelengths Returns: @@ -1261,69 +1539,94 @@ cpdef np.ndarray spring_assemble( for i in range(0, total_m): for j in range(0, total_m): # Submatrices for the initial stiffness - k11 = k_local[8 * i:8*i + 2, 8 * j:8*j + 2] - k12 = k_local[8 * i:8*i + 2, 8*j + 2:8*j + 4] - k13 = k_local[8 * i:8*i + 2, 8*j + 4:8*j + 6] - k14 = k_local[8 * i:8*i + 2, 8*j + 6:8*j + 8] - k21 = k_local[8*i + 2:8*i + 4, 8 * j:8*j + 2] - k22 = k_local[8*i + 2:8*i + 4, 8*j + 2:8*j + 4] - k23 = k_local[8*i + 2:8*i + 4, 8*j + 4:8*j + 6] - k24 = k_local[8*i + 2:8*i + 4, 8*j + 6:8*j + 8] - k31 = k_local[8*i + 4:8*i + 6, 8 * j:8*j + 2] - k32 = k_local[8*i + 4:8*i + 6, 8*j + 2:8*j + 4] - k33 = k_local[8*i + 4:8*i + 6, 8*j + 4:8*j + 6] - k34 = k_local[8*i + 4:8*i + 6, 8*j + 6:8*j + 8] - k41 = k_local[8*i + 6:8*i + 8, 8 * j:8*j + 2] - k42 = k_local[8*i + 6:8*i + 8, 8*j + 2:8*j + 4] - k43 = k_local[8*i + 6:8*i + 8, 8*j + 4:8*j + 6] - k44 = k_local[8*i + 6:8*i + 8, 8*j + 6:8*j + 8] - - k_global[4*n_nodes*i + (node_i+1) * 2 - 1:4*n_nodes*i + (node_i+1) * 2, - 4*n_nodes*j + (node_i+1) * 2 - 1:4*n_nodes*j + (node_i+1) * 2] += k11 + k11 = k_local[8 * i : 8 * i + 2, 8 * j : 8 * j + 2] + k12 = k_local[8 * i : 8 * i + 2, 8 * j + 2 : 8 * j + 4] + k13 = k_local[8 * i : 8 * i + 2, 8 * j + 4 : 8 * j + 6] + k14 = k_local[8 * i : 8 * i + 2, 8 * j + 6 : 8 * j + 8] + k21 = k_local[8 * i + 2 : 8 * i + 4, 8 * j : 8 * j + 2] + k22 = k_local[8 * i + 2 : 8 * i + 4, 8 * j + 2 : 8 * j + 4] + k23 = k_local[8 * i + 2 : 8 * i + 4, 8 * j + 4 : 8 * j + 6] + k24 = k_local[8 * i + 2 : 8 * i + 4, 8 * j + 6 : 8 * j + 8] + k31 = k_local[8 * i + 4 : 8 * i + 6, 8 * j : 8 * j + 2] + k32 = k_local[8 * i + 4 : 8 * i + 6, 8 * j + 2 : 8 * j + 4] + k33 = k_local[8 * i + 4 : 8 * i + 6, 8 * j + 4 : 8 * j + 6] + k34 = k_local[8 * i + 4 : 8 * i + 6, 8 * j + 6 : 8 * j + 8] + k41 = k_local[8 * i + 6 : 8 * i + 8, 8 * j : 8 * j + 2] + k42 = k_local[8 * i + 6 : 8 * i + 8, 8 * j + 2 : 8 * j + 4] + k43 = k_local[8 * i + 6 : 8 * i + 8, 8 * j + 4 : 8 * j + 6] + k44 = k_local[8 * i + 6 : 8 * i + 8, 8 * j + 6 : 8 * j + 8] + + k_global[ + 4 * n_nodes * i + (node_i + 1) * 2 - 1 : 4 * n_nodes * i + (node_i + 1) * 2, + 4 * n_nodes * j + (node_i + 1) * 2 - 1 : 4 * n_nodes * j + (node_i + 1) * 2, + ] += k11 if node_j != -1: - k_global[4*n_nodes*i + (node_i+1) * 2 - 1:4*n_nodes*i + (node_i+1) * 2, - 4*n_nodes*j + (node_j+1) * 2 - 1:4*n_nodes*j + (node_j+1) * 2] += k12 - k_global[4*n_nodes*i + (node_j+1) * 2 - 1:4*n_nodes*i + (node_j+1) * 2, - 4*n_nodes*j + (node_i+1) * 2 - 1:4*n_nodes*j + (node_i+1) * 2] += k21 - k_global[4*n_nodes*i + (node_j+1) * 2 - 1:4*n_nodes*i + (node_j+1) * 2, - 4*n_nodes*j + (node_j+1) * 2 - 1:4*n_nodes*j + (node_j+1) * 2] += k22 - - k_global[4*n_nodes*i + skip + (node_i+1) * 2 - 1:4*n_nodes*i + skip + (node_i+1) * 2, - 4*n_nodes*j + skip + (node_i+1) * 2 - 1:4*n_nodes*j + skip - + (node_i+1) * 2] += k33 + k_global[ + 4 * n_nodes * i + (node_i + 1) * 2 - 1 : 4 * n_nodes * i + (node_i + 1) * 2, + 4 * n_nodes * j + (node_j + 1) * 2 - 1 : 4 * n_nodes * j + (node_j + 1) * 2, + ] += k12 + k_global[ + 4 * n_nodes * i + (node_j + 1) * 2 - 1 : 4 * n_nodes * i + (node_j + 1) * 2, + 4 * n_nodes * j + (node_i + 1) * 2 - 1 : 4 * n_nodes * j + (node_i + 1) * 2, + ] += k21 + k_global[ + 4 * n_nodes * i + (node_j + 1) * 2 - 1 : 4 * n_nodes * i + (node_j + 1) * 2, + 4 * n_nodes * j + (node_j + 1) * 2 - 1 : 4 * n_nodes * j + (node_j + 1) * 2, + ] += k22 + + k_global[ + 4 * n_nodes * i + skip + (node_i + 1) * 2 - 1 : 4 * n_nodes * i + skip + (node_i + 1) * 2, + 4 * n_nodes * j + skip + (node_i + 1) * 2 - 1 : 4 * n_nodes * j + skip + (node_i + 1) * 2, + ] += k33 if node_j != -1: - k_global[4*n_nodes*i + skip + (node_i+1) * 2 - 1:4*n_nodes*i + skip - + (node_i+1) * 2, 4*n_nodes*j + skip + (node_j+1) * 2 - 1:4*n_nodes*j - + skip + (node_j+1) * 2] += k34 - k_global[4*n_nodes*i + skip + (node_j+1) * 2 - 1:4*n_nodes*i + skip - + (node_j+1) * 2, 4*n_nodes*j + skip + (node_i+1) * 2 - 1:4*n_nodes*j - + skip + (node_i+1) * 2] += k43 - k_global[4*n_nodes*i + skip + (node_j+1) * 2 - 1:4*n_nodes*i + skip - + (node_j+1) * 2, 4*n_nodes*j + skip + (node_j+1) * 2 - 1:4*n_nodes*j - + skip + (node_j+1) * 2] += k44 - - k_global[4*n_nodes*i + (node_i+1) * 2 - 1:4*n_nodes*i + (node_i+1) * 2, 4*n_nodes*j - + skip + (node_i+1) * 2 - 1:4*n_nodes*j + skip + (node_i+1) * 2] += k13 + k_global[ + 4 * n_nodes * i + skip + (node_i + 1) * 2 - 1 : 4 * n_nodes * i + skip + (node_i + 1) * 2, + 4 * n_nodes * j + skip + (node_j + 1) * 2 - 1 : 4 * n_nodes * j + skip + (node_j + 1) * 2, + ] += k34 + k_global[ + 4 * n_nodes * i + skip + (node_j + 1) * 2 - 1 : 4 * n_nodes * i + skip + (node_j + 1) * 2, + 4 * n_nodes * j + skip + (node_i + 1) * 2 - 1 : 4 * n_nodes * j + skip + (node_i + 1) * 2, + ] += k43 + k_global[ + 4 * n_nodes * i + skip + (node_j + 1) * 2 - 1 : 4 * n_nodes * i + skip + (node_j + 1) * 2, + 4 * n_nodes * j + skip + (node_j + 1) * 2 - 1 : 4 * n_nodes * j + skip + (node_j + 1) * 2, + ] += k44 + + k_global[ + 4 * n_nodes * i + (node_i + 1) * 2 - 1 : 4 * n_nodes * i + (node_i + 1) * 2, + 4 * n_nodes * j + skip + (node_i + 1) * 2 - 1 : 4 * n_nodes * j + skip + (node_i + 1) * 2, + ] += k13 if node_j != -1: - k_global[4*n_nodes*i + (node_i+1) * 2 - 1:4*n_nodes*i + (node_i+1) * 2, 4*n_nodes*j - + skip + (node_j+1) * 2 - 1:4*n_nodes*j + skip + (node_j+1) * 2] += k14 - k_global[4*n_nodes*i + (node_j+1) * 2 - 1:4*n_nodes*i + (node_j+1) * 2, 4*n_nodes*j - + skip + (node_i+1) * 2 - 1:4*n_nodes*j + skip + (node_i+1) * 2] += k23 - k_global[4*n_nodes*i + (node_j+1) * 2 - 1:4*n_nodes*i + (node_j+1) * 2, 4*n_nodes*j - + skip + (node_j+1) * 2 - 1:4*n_nodes*j + skip + (node_j+1) * 2] += k24 - - k_global[4*n_nodes*i + skip + (node_i+1) * 2 - 1:4*n_nodes*i + skip + (node_i+1) * 2, - 4*n_nodes*j + (node_i+1) * 2 - 1:4*n_nodes*j + (node_i+1) * 2] += k31 + k_global[ + 4 * n_nodes * i + (node_i + 1) * 2 - 1 : 4 * n_nodes * i + (node_i + 1) * 2, + 4 * n_nodes * j + skip + (node_j + 1) * 2 - 1 : 4 * n_nodes * j + skip + (node_j + 1) * 2, + ] += k14 + k_global[ + 4 * n_nodes * i + (node_j + 1) * 2 - 1 : 4 * n_nodes * i + (node_j + 1) * 2, + 4 * n_nodes * j + skip + (node_i + 1) * 2 - 1 : 4 * n_nodes * j + skip + (node_i + 1) * 2, + ] += k23 + k_global[ + 4 * n_nodes * i + (node_j + 1) * 2 - 1 : 4 * n_nodes * i + (node_j + 1) * 2, + 4 * n_nodes * j + skip + (node_j + 1) * 2 - 1 : 4 * n_nodes * j + skip + (node_j + 1) * 2, + ] += k24 + + k_global[ + 4 * n_nodes * i + skip + (node_i + 1) * 2 - 1 : 4 * n_nodes * i + skip + (node_i + 1) * 2, + 4 * n_nodes * j + (node_i + 1) * 2 - 1 : 4 * n_nodes * j + (node_i + 1) * 2, + ] += k31 if node_j != -1: - k_global[4*n_nodes*i + skip + (node_i+1) * 2 - 1:4*n_nodes*i + skip - + (node_i+1) * 2, - 4*n_nodes*j + (node_j+1) * 2 - 1:4*n_nodes*j + (node_j+1) * 2] += k32 - k_global[4*n_nodes*i + skip + (node_j+1) * 2 - 1:4*n_nodes*i + skip - + (node_j+1) * 2, - 4*n_nodes*j + (node_i+1) * 2 - 1:4*n_nodes*j + (node_i+1) * 2] += k41 - k_global[4*n_nodes*i + skip + (node_j+1) * 2 - 1:4*n_nodes*i + skip - + (node_j+1) * 2, - 4*n_nodes*j + (node_j+1) * 2 - 1:4*n_nodes*j + (node_j+1) * 2] += k42 + k_global[ + 4 * n_nodes * i + skip + (node_i + 1) * 2 - 1 : 4 * n_nodes * i + skip + (node_i + 1) * 2, + 4 * n_nodes * j + (node_j + 1) * 2 - 1 : 4 * n_nodes * j + (node_j + 1) * 2, + ] += k32 + k_global[ + 4 * n_nodes * i + skip + (node_j + 1) * 2 - 1 : 4 * n_nodes * i + skip + (node_j + 1) * 2, + 4 * n_nodes * j + (node_i + 1) * 2 - 1 : 4 * n_nodes * j + (node_i + 1) * 2, + ] += k41 + k_global[ + 4 * n_nodes * i + skip + (node_j + 1) * 2 - 1 : 4 * n_nodes * i + skip + (node_j + 1) * 2, + 4 * n_nodes * j + (node_j + 1) * 2 - 1 : 4 * n_nodes * j + (node_j + 1) * 2, + ] += k42 return k_global @@ -1342,21 +1645,20 @@ cpdef double ym_at_ys(str b_c, double m_i, double y_s, double length): Returns: y_m (float): longitudinal shape function value - + BWS in 2015 """ cdef double y_m - if b_c == 'S-S': + if b_c == "S-S": y_m = np.sin(m_i * np.pi * y_s / length) - elif b_c == 'C-C': + elif b_c == "C-C": y_m = np.sin(m_i * np.pi * y_s / length) * np.sin(np.pi * y_s / length) - elif b_c == 'S-C' or b_c == 'C-S': - y_m = np.sin((m_i+1) * np.pi * y_s / length - ) + (m_i+1) / m_i * np.sin(m_i * np.pi * y_s / length) - elif b_c == 'C-F' or b_c == 'F-C': - y_m = 1 - np.cos((m_i-0.5) * np.pi * y_s / length) - elif b_c == 'C-G' or b_c == 'G-C': - y_m = np.sin((m_i-0.5) * np.pi * y_s / length) * np.sin(np.pi * y_s / length / 2) + elif b_c in ("S-C", "C-S"): + y_m = np.sin((m_i + 1) * np.pi * y_s / length) + (m_i + 1) / m_i * np.sin(m_i * np.pi * y_s / length) + elif b_c in ("C-F", "F-C"): + y_m = 1 - np.cos((m_i - 0.5) * np.pi * y_s / length) + elif b_c in ("C-G", "G-C"): + y_m = np.sin((m_i - 0.5) * np.pi * y_s / length) * np.sin(np.pi * y_s / length / 2) else: raise ValueError(f"Unrecognised boundary condition '{b_c}'") @@ -1381,23 +1683,27 @@ cpdef double ymprime_at_ys(str b_c, double m_i, double y_s, double length): BWS in 2015 """ cdef double y_m_prime - if b_c == 'S-S': + if b_c == "S-S": y_m_prime = (np.pi * m_i * np.cos((np.pi * m_i * y_s) / length)) / length - elif b_c == 'C-C': - y_m_prime = (np.pi * np.cos((np.pi*y_s) / length) \ - * np.sin((np.pi*m_i*y_s) / length)) / length \ - + (np.pi*m_i * np.sin((np.pi*y_s)/length) \ - * np.cos((np.pi*m_i*y_s)/length)) / length - elif b_c == 'S-C' or b_c == 'C-S': - y_m_prime = (np.pi * np.cos((np.pi*y_s * (m_i + 1))/length) * (m_i + 1)) / length \ - + (np.pi * np.cos((np.pi*m_i*y_s)/length)*(m_i + 1)) / length - elif b_c == 'C-F' or b_c == 'F-C': - y_m_prime = (np.pi * np.sin((np.pi * y_s * (m_i - 1/2)) / length) * (m_i - 1/2)) / length - elif b_c == 'C-G' or b_c == 'G-C': - y_m_prime = (np.pi*np.sin((np.pi*y_s * (m_i - 1/2))/length) \ - * np.cos((np.pi*y_s)/(2*length)))/(2*length) \ - + (np.pi*np.cos((np.pi*y_s*(m_i - 1/2))/length) \ - * np.sin((np.pi*y_s)/(2*length))*(m_i - 1/2))/length + elif b_c == "C-C": + y_m_prime = (np.pi * np.cos((np.pi * y_s) / length) * np.sin((np.pi * m_i * y_s) / length)) / length + ( + np.pi * m_i * np.sin((np.pi * y_s) / length) * np.cos((np.pi * m_i * y_s) / length) + ) / length + elif b_c in ("S-C", "C-S"): + y_m_prime = (np.pi * np.cos((np.pi * y_s * (m_i + 1)) / length) * (m_i + 1)) / length + ( + np.pi * np.cos((np.pi * m_i * y_s) / length) * (m_i + 1) + ) / length + elif b_c in ("C-F", "F-C"): + y_m_prime = (np.pi * np.sin((np.pi * y_s * (m_i - 1 / 2)) / length) * (m_i - 1 / 2)) / length + elif b_c in ("C-G", "G-C"): + y_m_prime = (np.pi * np.sin((np.pi * y_s * (m_i - 1 / 2)) / length) * np.cos((np.pi * y_s) / (2 * length))) / ( + 2 * length + ) + ( + np.pi + * np.cos((np.pi * y_s * (m_i - 1 / 2)) / length) + * np.sin((np.pi * y_s) / (2 * length)) + * (m_i - 1 / 2) + ) / length else: raise ValueError(f"Unrecognised boundary condition '{b_c}'") diff --git a/pycufsm/analysis_p.py b/pycufsm/analysis_p.py index 121633e..d49a032 100644 --- a/pycufsm/analysis_p.py +++ b/pycufsm/analysis_p.py @@ -256,11 +256,27 @@ def k_kg_local( i_5=i_5, ) k_local[8 * i + 4 : 8 * (i + 1), 8 * j + 4 : 8 * (j + 1)] = calc_kf_mp( - d_x=d_x, d_y=d_y, d_1=d_1, d_xy=d_xy, b_strip=b_strip, i_1=i_1, i_2=i_2, i_3=i_3, i_4=i_4, i_5=i_5 + d_x=d_x, + d_y=d_y, + d_1=d_1, + d_xy=d_xy, + b_strip=b_strip, + i_1=i_1, + i_2=i_2, + i_3=i_3, + i_4=i_4, + i_5=i_5, ) kg_local[8 * i : 8 * i + 4, 8 * j : 8 * j + 4] = calc_gm_mp( - u_i=u_i, u_j=u_j, b_strip=b_strip, length=length, ty_1=ty_1, ty_2=ty_2, i_4=i_4, i_5=i_5 + u_i=u_i, + u_j=u_j, + b_strip=b_strip, + length=length, + ty_1=ty_1, + ty_2=ty_2, + i_4=i_4, + i_5=i_5, ) kg_local[8 * i + 4 : 8 * (i + 1), 8 * j + 4 : 8 * (j + 1)] = calc_gf_mp( ty_1=ty_1, ty_2=ty_2, b_strip=b_strip, i_5=i_5 @@ -338,7 +354,11 @@ def kglobal_transv( node_i = int(elements[i, 1]) node_j = int(elements[i, 2]) k_global_transv = assemble_single( - k_global=k_global_transv, k_local=k_local, node_i=node_i, node_j=node_j, n_nodes=n_nodes + k_global=k_global_transv, + k_local=k_local, + node_i=node_i, + node_j=node_j, + n_nodes=n_nodes, ) return k_global_transv @@ -396,10 +416,10 @@ def klocal_transv( c_2 = u_p / length [i_1, _, _, _, _] = bc_i1_5(b_c=b_c, m_i=m_i, m_j=m_i, length=length) - i_2 = 0 - i_3 = 0 - i_4 = 0 - i_5 = 0 + i_2 = 0.0 + i_3 = 0.0 + i_4 = 0.0 + i_5 = 0.0 k_local[0:4, 0:4] = calc_km_mp( e_1=e_1, @@ -417,7 +437,16 @@ def klocal_transv( i_5=i_5, ) k_local[4:8, 4:8] = calc_kf_mp( - d_x=d_x, d_y=d_y, d_1=d_1, d_xy=d_xy, b_strip=b_strip, i_1=i_1, i_2=i_2, i_3=i_3, i_4=i_4, i_5=i_5 + d_x=d_x, + d_y=d_y, + d_1=d_1, + d_xy=d_xy, + b_strip=b_strip, + i_1=i_1, + i_2=i_2, + i_3=i_3, + i_4=i_4, + i_5=i_5, ) return k_local @@ -671,7 +700,14 @@ def calc_kf_mp( def calc_gm_mp( - u_i: float, u_j: float, b_strip: float, length: float, ty_1: float, ty_2: float, i_4: float, i_5: float + u_i: float, + u_j: float, + b_strip: float, + length: float, + ty_1: float, + ty_2: float, + i_4: float, + i_5: float, ) -> np.ndarray: """Calculate the membrane geometric stiffness sub-matrix, used in the assembly of local geometric stiffness matrices @@ -750,7 +786,7 @@ def calc_gf_mp(ty_1: float, ty_2: float, b_strip: float, i_5: float) -> np.ndarr return gf_mp -def bc_i1_5(b_c: str, m_i: float, m_j: float, length: float) -> list: +def bc_i1_5(b_c: str, m_i: float, m_j: float, length: float) -> Tuple[float, float, float, float, float]: """Calculate the 5 undetermined parameters i_1,i_2,i_3,i_4,i_5 for local elastic and geometric stiffness matrices. @@ -884,7 +920,7 @@ def bc_i1_5(b_c: str, m_i: float, m_j: float, length: float) -> list: i_4 = -(m_i**4) * np.pi**4 / 8 / length**3 i_5 = -(m_i**2) * np.pi**2 / 8 / length - return [i_1, i_2, i_3, i_4, i_5] + return (i_1, i_2, i_3, i_4, i_5) def trans(alpha: float, total_m: int) -> np.ndarray: @@ -1137,7 +1173,13 @@ def assemble( return k_global, kg_global -def assemble_single(k_global: np.ndarray, k_local: np.ndarray, node_i: int, node_j: int, n_nodes: int) -> np.ndarray: +def assemble_single( + k_global: np.ndarray, + k_local: np.ndarray, + node_i: int, + node_j: int, + n_nodes: int, +) -> np.ndarray: """this routine adds the element contribution to the global stiffness matrix basically it does the same as routine 'assemble', however: it does not care about kg_global (geom stiff matrix) @@ -1181,10 +1223,22 @@ def assemble_single(k_global: np.ndarray, k_local: np.ndarray, node_i: int, node k_global[node_j * 2 : node_j * 2 + 2, node_i * 2 : node_i * 2 + 2] += k21 k_global[node_j * 2 : node_j * 2 + 2, node_j * 2 : node_j * 2 + 2] += k22 - k_global[skip + node_i * 2 : skip + node_i * 2 + 2, skip + node_i * 2 : skip + node_i * 2 + 2] += k33 - k_global[skip + node_i * 2 : skip + node_i * 2 + 2, skip + node_j * 2 : skip + node_j * 2 + 2] += k34 - k_global[skip + node_j * 2 : skip + node_j * 2 + 2, skip + node_i * 2 : skip + node_i * 2 + 2] += k43 - k_global[skip + node_j * 2 : skip + node_j * 2 + 2, skip + node_j * 2 : skip + node_j * 2 + 2] += k44 + k_global[ + skip + node_i * 2 : skip + node_i * 2 + 2, + skip + node_i * 2 : skip + node_i * 2 + 2, + ] += k33 + k_global[ + skip + node_i * 2 : skip + node_i * 2 + 2, + skip + node_j * 2 : skip + node_j * 2 + 2, + ] += k34 + k_global[ + skip + node_j * 2 : skip + node_j * 2 + 2, + skip + node_i * 2 : skip + node_i * 2 + 2, + ] += k43 + k_global[ + skip + node_j * 2 : skip + node_j * 2 + 2, + skip + node_j * 2 : skip + node_j * 2 + 2, + ] += k44 k_global[node_i * 2 : node_i * 2 + 2, skip + node_i * 2 : skip + node_i * 2 + 2] += k13 k_global[node_i * 2 : node_i * 2 + 2, skip + node_j * 2 : skip + node_j * 2 + 2] += k14 @@ -1200,7 +1254,15 @@ def assemble_single(k_global: np.ndarray, k_local: np.ndarray, node_i: int, node def spring_klocal( - k_u: float, k_v: float, k_w: float, k_q: float, length: float, b_c: str, m_a: np.ndarray, discrete: int, y_s: float + k_u: float, + k_v: float, + k_w: float, + k_q: float, + length: float, + b_c: str, + m_a: np.ndarray, + discrete: int, + y_s: float, ) -> np.ndarray: """Generate spring stiffness matrix (k_local) in local coordinates, modified from klocal @@ -1245,9 +1307,19 @@ def spring_klocal( km_mp = np.array( [ [k_u * i_1, 0, -k_u * i_1, 0], - [0, k_v * i_5 * length**2 / (u_i * u_j), 0, -k_v * i_5 * length**2 / (u_i * u_j)], + [ + 0, + k_v * i_5 * length**2 / (u_i * u_j), + 0, + -k_v * i_5 * length**2 / (u_i * u_j), + ], [-k_u * i_1, 0, k_u * i_1, 0], - [0, -k_v * i_5 * length**2 / (u_i * u_j), 0, k_v * i_5 * length**2 / (u_i * u_j)], + [ + 0, + -k_v * i_5 * length**2 / (u_i * u_j), + 0, + k_v * i_5 * length**2 / (u_i * u_j), + ], ] ) # assemble the matrix of kf_mp (flexural stiffness) @@ -1296,7 +1368,12 @@ def bc_i1_5_atpoint(b_c: str, m_i: float, m_j: float, length: float, y_s: float) def spring_assemble( - k_global: np.ndarray, k_local: np.ndarray, node_i: int, node_j: int, n_nodes: int, m_a: np.ndarray + k_global: np.ndarray, + k_local: np.ndarray, + node_i: int, + node_j: int, + n_nodes: int, + m_a: np.ndarray, ) -> np.ndarray: """Add the (spring) contribution to the global stiffness matrix diff --git a/pycufsm/cfsm.py b/pycufsm/cfsm.py index 3aa79ba..655ba63 100644 --- a/pycufsm/cfsm.py +++ b/pycufsm/cfsm.py @@ -197,9 +197,9 @@ def base_update( for i, m_i in enumerate(m_a): b_v_m = b_v_l[n_dof_m * i : n_dof_m * (i + 1), n_dof_m * i : n_dof_m * (i + 1)] # k_global/kg_global - if gbt_con["norm"] in (2, 3) or gbt_con["o_space"] in (2, 3) or gbt_con["orth"] in (2, 3): + if gbt_con["norm"] in (2, 3) or gbt_con["o_space"] in (2, 3, 4) or gbt_con["orth"] in (2, 3): # axial loading or real loading by either gbt_con['orth'] = 2 or gbt_con['orth'] = 3 - if gbt_con["orth"] == 1 or gbt_con["orth"] == 2: + if gbt_con["orth"] in (1, 2): nodes_base = deepcopy(nodes) nodes_base[:, 7] = np.ones_like(nodes[:, 7]) # set u_p stress to 1.0 (axial) k_global, kg_global = analysis.k_kg_global( @@ -224,13 +224,7 @@ def base_update( # orthogonalization/normalization begins # - if ( - gbt_con["orth"] == 2 - or gbt_con["orth"] == 3 - or gbt_con["o_space"] == 2 - or gbt_con["o_space"] == 3 - or gbt_con["o_space"] == 4 - ): + if gbt_con["orth"] in (2, 3) or gbt_con["o_space"] in (2, 3, 4): # indices if gbt_con["o_space"] == 1: dof_index = np.zeros((5, 2)) @@ -267,19 +261,18 @@ def base_update( if dof_sub[1] >= dof_sub[0]: k_global_sub = b_v_m[:, dof_sub0:dof_sub1].conj().T @ k_global @ b_v_m[:, dof_sub0:dof_sub1] kg_global_sub = b_v_m[:, dof_sub0:dof_sub1].conj().T @ kg_global @ b_v_m[:, dof_sub0:dof_sub1] - [eigenvalues, eigenvectors] = spla.eig(a=k_global_sub, b=kg_global_sub) + eigenvalues, eigenvectors = spla.eig(a=k_global_sub, b=kg_global_sub) lf_sub = np.real(eigenvalues) indexsub = np.argsort(lf_sub) lf_sub = lf_sub[indexsub] eigenvectors = np.real(eigenvectors[:, indexsub]) - if gbt_con["norm"] == 2 or gbt_con["norm"] == 3: + if gbt_con["norm"] in (2, 3): if gbt_con["norm"] == 2: - s_matrix = eigenvectors.conj().T @ k_global_sub @ eigenvectors + s_matrix = np.diag(eigenvectors.conj().T @ k_global_sub @ eigenvectors) - if gbt_con["norm"] == 3: - s_matrix = eigenvectors.conj().T @ kg_global_sub @ eigenvectors + elif gbt_con["norm"] == 3: + s_matrix = np.diag(eigenvectors.conj().T @ kg_global_sub @ eigenvectors) - s_matrix = np.diag(s_matrix) for j in range(0, int(dof_sub[1] - dof_sub[0])): eigenvectors[:, j] = np.transpose( np.conj(np.linalg.lstsq(eigenvectors[:, j].conj().T, np.sqrt(s_matrix).conj().T)) @@ -288,7 +281,7 @@ def base_update( b_v_m[:, dof_sub0:dof_sub1] = b_v_m[:, dof_sub0:dof_sub1] @ eigenvectors # normalization for gbt_con['o_space'] = 1 - if (gbt_con["norm"] == 2 or gbt_con["norm"] == 3) and gbt_con["o_space"] == 1: + if gbt_con["norm"] in (2, 3) and gbt_con["o_space"] == 1: for j in range(0, n_dof_m): if gbt_con["norm"] == 2: b_v_m[:, j] = np.transpose( @@ -317,12 +310,12 @@ def base_update( b_v[n_dof_m * i : n_dof_m * (i + 1), n_dof_m * i : n_dof_m * (i + 1)] = b_v_m - else: + elif gbt_con["couple"] != 1: # coupled basis # k_global/kg_global if gbt_con["norm"] in (2, 3) or gbt_con["o_space"] in (2, 3) or gbt_con["orth"] in (2, 3): # axial loading or real loading by either gbt_con['orth'] = 2 or gbt_con['orth'] = 3 - if gbt_con["orth"] == 1 or gbt_con["orth"] == 2: + if gbt_con["orth"] in (1, 2): nodes_base = deepcopy(nodes) nodes_base[:, 7] = np.ones_like(nodes[:, 7]) # set u_p stress to 1.0 (axial) else: @@ -333,13 +326,7 @@ def base_update( ) # orthogonalization/normalization begins - if ( - gbt_con["orth"] == 2 - or gbt_con["orth"] == 3 - or gbt_con["o_space"] == 2 - or gbt_con["o_space"] == 3 - or gbt_con["o_space"] == 4 - ): + if gbt_con["orth"] in (2, 3) or gbt_con["o_space"] in (2, 3, 4): # indices dof_index = np.zeros((4, 2)) dof_index[0, 0] = 0 @@ -420,7 +407,7 @@ def base_update( indexsub = np.argsort(lf_sub) lf_sub = lf_sub[indexsub] eigenvectors = np.real(eigenvectors[:, indexsub]) - if gbt_con["norm"] == 2 or gbt_con["norm"] == 3: + if gbt_con["norm"] in (2, 3): if gbt_con["norm"] == 2: s_matrix = eigenvectors.conj().T @ k_global_sub @ eigenvectors if gbt_con["norm"] == 3: @@ -459,7 +446,7 @@ def base_update( ] # normalization for gbt_con['o_space'] = 1 - if (gbt_con["norm"] == 2 or gbt_con["norm"] == 3) and (gbt_con["o_space"] == 1): + if gbt_con["norm"] in (2, 3) and gbt_con["o_space"] == 1: for i in range(0, n_dof_m * total_m): if gbt_con["norm"] == 2: b_v[:, i] = np.transpose( @@ -1834,7 +1821,7 @@ def dof_ordering(node_props: np.ndarray) -> np.ndarray: i_c = 0 i_s = 0 for i, n_prop in enumerate(node_props): - if n_prop[3] == 1 or n_prop[3] == 2: # corner or edge nodes + if n_prop[3] in (1, 2): # corner or edge nodes dof_perm[2 * i + 1, i_c] = 1 i_c = i_c + 1 if n_prop[3] == 3: # sub nodes @@ -1860,7 +1847,7 @@ def dof_ordering(node_props: np.ndarray) -> np.ndarray: i_c = 0 i_s = 0 for i, n_prop in enumerate(node_props): - if n_prop[3] == 1 or n_prop[3] == 2: # corner or edge nodes + if n_prop[3] in (1, 2): # corner or edge nodes dof_perm[2 * n_node_props + 2 * i + 1, 3 * n_main_nodes + i_c] = 1 i_c = i_c + 1 if n_prop[3] == 3: # sub nodes